Eigen-Contrib  5.0.1
 
Loading...
Searching...
No Matches
fftw_impl.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009 Mark Borgerding mark a borgerding net
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9// SPDX-License-Identifier: MPL-2.0
10
11#ifndef EIGEN_FFT_FFTW_IMPL_H
12#define EIGEN_FFT_FFTW_IMPL_H
13
14// IWYU pragma: private
15#include "./InternalHeaderCheck.h"
16
17#include <memory>
18#include <mutex>
19
20namespace Eigen {
21
22namespace internal {
23
24// FFTW uses non-const arguments,
25// so const_cast is needed for all the args it uses.
26//
27// This should be safe as long as
28// 1. we use FFTW_ESTIMATE for all our planning
29// see the FFTW docs section 4.3.2 "Planner Flags"
30// 2. fftw_complex is compatible with std::complex
31// This assumes std::complex<T> layout is array of size 2 with real,imag
32template <typename T>
33inline T *fftw_cast(const T *p) {
34 return const_cast<T *>(p);
35}
36
37inline fftw_complex *fftw_cast(const std::complex<double> *p) {
38 return const_cast<fftw_complex *>(reinterpret_cast<const fftw_complex *>(p));
39}
40
41inline fftwf_complex *fftw_cast(const std::complex<float> *p) {
42 return const_cast<fftwf_complex *>(reinterpret_cast<const fftwf_complex *>(p));
43}
44
45inline fftwl_complex *fftw_cast(const std::complex<long double> *p) {
46 return const_cast<fftwl_complex *>(reinterpret_cast<const fftwl_complex *>(p));
47}
48
49// The FFTW planner is not thread-safe: fftw_execute and its new-array variants,
50// which is what this backend runs transforms through, are the only entry points
51// that may be called concurrently (FFTW manual, "Thread safety"), so plan
52// creation and destruction are serialized through this mutex. A template
53// static gives the header-only definition; std::mutex is
54// constexpr-constructible, so the mutex is ready before any thread starts.
55// The planner state it stands in for is one per process, so the mutex must be
56// too: the explicit default visibility is what keeps the definition from being
57// bound locally under -fvisibility=hidden, where each library planning through
58// Eigen would get a mutex of its own and serialize nothing between them. The
59// module documentation covers what no header can reach, which needs FFTW's own
60// fftw_make_planner_thread_safe().
61#if EIGEN_HAS_ATTRIBUTE(visibility) && !EIGEN_OS_WIN
62#define EIGEN_FFTW_PLANNER_MUTEX_VISIBILITY __attribute__((visibility("default")))
63#else
64#define EIGEN_FFTW_PLANNER_MUTEX_VISIBILITY
65#endif
66
67template <typename Dummy = void>
68struct fftw_planner_lock {
69 static EIGEN_FFTW_PLANNER_MUTEX_VISIBILITY std::mutex mutex;
70};
71template <typename Dummy>
72EIGEN_FFTW_PLANNER_MUTEX_VISIBILITY std::mutex fftw_planner_lock<Dummy>::mutex;
73
74inline std::mutex &fftw_planner_mutex() { return fftw_planner_lock<>::mutex; }
75
76template <typename PlanFactory>
77inline decltype(auto) fftw_make_plan(PlanFactory factory) {
78 std::lock_guard<std::mutex> lock(fftw_planner_mutex());
79 return factory();
80}
81
82template <typename T>
83struct fftw_plan {};
84
85template <>
86struct fftw_plan<float> {
87 typedef float scalar_type;
88 typedef fftwf_complex complex_type;
89 std::shared_ptr<fftwf_plan_s> m_plan;
90 fftw_plan() = default;
91
92 void set_plan(fftwf_plan p) {
93 m_plan.reset(p, [](fftwf_plan plan) {
94 std::lock_guard<std::mutex> lock(fftw_planner_mutex());
95 fftwf_destroy_plan(plan);
96 });
97 }
98 inline void fwd(complex_type *dst, complex_type *src, int nfft) {
99 if (!m_plan)
100 set_plan(fftw_make_plan(
101 [&] { return fftwf_plan_dft_1d(nfft, src, dst, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
102 fftwf_execute_dft(m_plan.get(), src, dst);
103 }
104 inline void inv(complex_type *dst, complex_type *src, int nfft) {
105 if (!m_plan)
106 set_plan(fftw_make_plan(
107 [&] { return fftwf_plan_dft_1d(nfft, src, dst, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
108 fftwf_execute_dft(m_plan.get(), src, dst);
109 }
110 inline void fwd(complex_type *dst, scalar_type *src, int nfft) {
111 if (!m_plan)
112 set_plan(
113 fftw_make_plan([&] { return fftwf_plan_dft_r2c_1d(nfft, src, dst, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
114 fftwf_execute_dft_r2c(m_plan.get(), src, dst);
115 }
116 inline void inv(scalar_type *dst, complex_type *src, int nfft) {
117 if (!m_plan)
118 set_plan(
119 fftw_make_plan([&] { return fftwf_plan_dft_c2r_1d(nfft, src, dst, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
120 fftwf_execute_dft_c2r(m_plan.get(), src, dst);
121 }
122
123 inline void fwd2(complex_type *dst, complex_type *src, int n0, int n1) {
124 if (!m_plan)
125 set_plan(fftw_make_plan(
126 [&] { return fftwf_plan_dft_2d(n0, n1, src, dst, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
127 fftwf_execute_dft(m_plan.get(), src, dst);
128 }
129 inline void inv2(complex_type *dst, complex_type *src, int n0, int n1) {
130 if (!m_plan)
131 set_plan(fftw_make_plan(
132 [&] { return fftwf_plan_dft_2d(n0, n1, src, dst, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
133 fftwf_execute_dft(m_plan.get(), src, dst);
134 }
135};
136template <>
137struct fftw_plan<double> {
138 typedef double scalar_type;
139 typedef fftw_complex complex_type;
140 std::shared_ptr<fftw_plan_s> m_plan;
141 fftw_plan() = default;
142
143 void set_plan(::fftw_plan p) {
144 m_plan.reset(p, [](::fftw_plan plan) {
145 std::lock_guard<std::mutex> lock(fftw_planner_mutex());
146 fftw_destroy_plan(plan);
147 });
148 }
149 inline void fwd(complex_type *dst, complex_type *src, int nfft) {
150 if (!m_plan)
151 set_plan(fftw_make_plan(
152 [&] { return fftw_plan_dft_1d(nfft, src, dst, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
153 fftw_execute_dft(m_plan.get(), src, dst);
154 }
155 inline void inv(complex_type *dst, complex_type *src, int nfft) {
156 if (!m_plan)
157 set_plan(fftw_make_plan(
158 [&] { return fftw_plan_dft_1d(nfft, src, dst, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
159 fftw_execute_dft(m_plan.get(), src, dst);
160 }
161 inline void fwd(complex_type *dst, scalar_type *src, int nfft) {
162 if (!m_plan)
163 set_plan(
164 fftw_make_plan([&] { return fftw_plan_dft_r2c_1d(nfft, src, dst, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
165 fftw_execute_dft_r2c(m_plan.get(), src, dst);
166 }
167 inline void inv(scalar_type *dst, complex_type *src, int nfft) {
168 if (!m_plan)
169 set_plan(
170 fftw_make_plan([&] { return fftw_plan_dft_c2r_1d(nfft, src, dst, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
171 fftw_execute_dft_c2r(m_plan.get(), src, dst);
172 }
173 inline void fwd2(complex_type *dst, complex_type *src, int n0, int n1) {
174 if (!m_plan)
175 set_plan(fftw_make_plan(
176 [&] { return fftw_plan_dft_2d(n0, n1, src, dst, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
177 fftw_execute_dft(m_plan.get(), src, dst);
178 }
179 inline void inv2(complex_type *dst, complex_type *src, int n0, int n1) {
180 if (!m_plan)
181 set_plan(fftw_make_plan(
182 [&] { return fftw_plan_dft_2d(n0, n1, src, dst, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
183 fftw_execute_dft(m_plan.get(), src, dst);
184 }
185};
186template <>
187struct fftw_plan<long double> {
188 typedef long double scalar_type;
189 typedef fftwl_complex complex_type;
190 std::shared_ptr<fftwl_plan_s> m_plan;
191 fftw_plan() = default;
192
193 void set_plan(fftwl_plan p) {
194 m_plan.reset(p, [](fftwl_plan plan) {
195 std::lock_guard<std::mutex> lock(fftw_planner_mutex());
196 fftwl_destroy_plan(plan);
197 });
198 }
199 inline void fwd(complex_type *dst, complex_type *src, int nfft) {
200 if (!m_plan)
201 set_plan(fftw_make_plan(
202 [&] { return fftwl_plan_dft_1d(nfft, src, dst, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
203 fftwl_execute_dft(m_plan.get(), src, dst);
204 }
205 inline void inv(complex_type *dst, complex_type *src, int nfft) {
206 if (!m_plan)
207 set_plan(fftw_make_plan(
208 [&] { return fftwl_plan_dft_1d(nfft, src, dst, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
209 fftwl_execute_dft(m_plan.get(), src, dst);
210 }
211 inline void fwd(complex_type *dst, scalar_type *src, int nfft) {
212 if (!m_plan)
213 set_plan(
214 fftw_make_plan([&] { return fftwl_plan_dft_r2c_1d(nfft, src, dst, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
215 fftwl_execute_dft_r2c(m_plan.get(), src, dst);
216 }
217 inline void inv(scalar_type *dst, complex_type *src, int nfft) {
218 if (!m_plan)
219 set_plan(
220 fftw_make_plan([&] { return fftwl_plan_dft_c2r_1d(nfft, src, dst, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
221 fftwl_execute_dft_c2r(m_plan.get(), src, dst);
222 }
223 inline void fwd2(complex_type *dst, complex_type *src, int n0, int n1) {
224 if (!m_plan)
225 set_plan(fftw_make_plan(
226 [&] { return fftwl_plan_dft_2d(n0, n1, src, dst, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
227 fftwl_execute_dft(m_plan.get(), src, dst);
228 }
229 inline void inv2(complex_type *dst, complex_type *src, int n0, int n1) {
230 if (!m_plan)
231 set_plan(fftw_make_plan(
232 [&] { return fftwl_plan_dft_2d(n0, n1, src, dst, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); }));
233 fftwl_execute_dft(m_plan.get(), src, dst);
234 }
235};
236
237template <typename Scalar_>
238struct fftw_impl {
239 typedef Scalar_ Scalar;
240 typedef std::complex<Scalar> Complex;
241
242 inline void clear() { m_plans.clear(); }
243
244 // complex-to-complex forward FFT
245 inline void fwd(Complex *dst, const Complex *src, int nfft) {
246 get_plan(nfft, false, /*real_io=*/false, dst, src).fwd(fftw_cast(dst), fftw_cast(src), nfft);
247 }
248
249 // real-to-complex forward FFT
250 inline void fwd(Complex *dst, const Scalar *src, int nfft) {
251 get_plan(nfft, false, /*real_io=*/true, dst, src).fwd(fftw_cast(dst), fftw_cast(src), nfft);
252 }
253
254 // 2-d complex-to-complex
255 inline void fwd2(Complex *dst, const Complex *src, int n0, int n1) {
256 get_plan(n0, n1, false, /*real_io=*/false, dst, src).fwd2(fftw_cast(dst), fftw_cast(src), n0, n1);
257 }
258
259 // inverse complex-to-complex
260 inline void inv(Complex *dst, const Complex *src, int nfft) {
261 get_plan(nfft, true, /*real_io=*/false, dst, src).inv(fftw_cast(dst), fftw_cast(src), nfft);
262 }
263
264 // half-complex to scalar
265 inline void inv(Scalar *dst, const Complex *src, int nfft) {
266 get_plan(nfft, true, /*real_io=*/true, dst, src).inv(fftw_cast(dst), fftw_cast(src), nfft);
267 }
268
269 // 2-d complex-to-complex
270 inline void inv2(Complex *dst, const Complex *src, int n0, int n1) {
271 get_plan(n0, n1, true, /*real_io=*/false, dst, src).inv2(fftw_cast(dst), fftw_cast(src), n0, n1);
272 }
273
274 protected:
275 typedef fftw_plan<Scalar> PlanData;
276
277 typedef Eigen::numext::int64_t int64_t;
278
279 typedef std::map<int64_t, PlanData> PlanMap;
280
281 PlanMap m_plans;
282
283 // Pack (inverse, real_io, inplace, aligned) into 4 contiguous low bits of
284 // the cache key. real_io distinguishes r2c/c2r from c2c so that reusing
285 // the same FFT object across real-input and complex-input transforms
286 // doesn't return a cached plan of the wrong kind.
287 static int64_t plan_flags(bool inverse, bool real_io, void *dst, const void *src) {
288 bool inplace = (dst == src);
289 bool aligned = ((reinterpret_cast<size_t>(src) & 15) | (reinterpret_cast<size_t>(dst) & 15)) == 0;
290 return (inverse << 3) | (real_io << 2) | (inplace << 1) | aligned;
291 }
292
293 inline PlanData &get_plan(int nfft, bool inverse, bool real_io, void *dst, const void *src) {
294 int64_t key = ((nfft << 4) | plan_flags(inverse, real_io, dst, src)) << 1;
295 return m_plans[key];
296 }
297
298 inline PlanData &get_plan(int n0, int n1, bool inverse, bool real_io, void *dst, const void *src) {
299 int64_t key = (((((int64_t)n0) << 31) | (n1 << 4) | plan_flags(inverse, real_io, dst, src)) << 1) + 1;
300 return m_plans[key];
301 }
302};
303
304} // end namespace internal
305
306} // end namespace Eigen
307
308#endif // EIGEN_FFT_FFTW_IMPL_H
Eigen
Namespace containing all symbols from the Eigen library.