ei_fftw_impl.h
Go to the documentation of this file.
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 
10 #include "./InternalHeaderCheck.h"
11 
12 #include <memory>
13 
14 namespace Eigen {
15 
16 namespace internal {
17 
18  // FFTW uses non-const arguments
19  // so we must use ugly const_cast calls for all the args it uses
20  //
21  // This should be safe as long as
22  // 1. we use FFTW_ESTIMATE for all our planning
23  // see the FFTW docs section 4.3.2 "Planner Flags"
24  // 2. fftw_complex is compatible with std::complex
25  // This assumes std::complex<T> layout is array of size 2 with real,imag
26  template <typename T>
27  inline
28  T * fftw_cast(const T* p)
29  {
30  return const_cast<T*>( p);
31  }
32 
33  inline
34  fftw_complex * fftw_cast( const std::complex<double> * p)
35  {
36  return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) );
37  }
38 
39  inline
40  fftwf_complex * fftw_cast( const std::complex<float> * p)
41  {
42  return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) );
43  }
44 
45  inline
46  fftwl_complex * fftw_cast( const std::complex<long double> * p)
47  {
48  return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) );
49  }
50 
51  template <typename T>
52  struct fftw_plan {};
53 
54  template <>
55  struct fftw_plan<float>
56  {
57  typedef float scalar_type;
58  typedef fftwf_complex complex_type;
59  std::shared_ptr<fftwf_plan_s> m_plan;
60  fftw_plan() = default;
61 
62  void set_plan(fftwf_plan p) { m_plan.reset(p, fftwf_destroy_plan); }
63  inline
64  void fwd(complex_type * dst,complex_type * src,int nfft) {
65  if (m_plan==NULL) set_plan(fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
66  fftwf_execute_dft( m_plan.get(), src,dst);
67  }
68  inline
69  void inv(complex_type * dst,complex_type * src,int nfft) {
70  if (m_plan==NULL) set_plan(fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
71  fftwf_execute_dft( m_plan.get(), src,dst);
72  }
73  inline
74  void fwd(complex_type * dst,scalar_type * src,int nfft) {
75  if (m_plan==NULL) set_plan(fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
76  fftwf_execute_dft_r2c( m_plan.get(),src,dst);
77  }
78  inline
79  void inv(scalar_type * dst,complex_type * src,int nfft) {
80  if (m_plan==NULL)
81  set_plan(fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
82  fftwf_execute_dft_c2r( m_plan.get(), src,dst);
83  }
84 
85  inline
86  void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
87  if (m_plan==NULL) set_plan(fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
88  fftwf_execute_dft( m_plan.get(), src,dst);
89  }
90  inline
91  void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
92  if (m_plan==NULL) set_plan(fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
93  fftwf_execute_dft( m_plan.get(), src,dst);
94  }
95 
96  };
97  template <>
98  struct fftw_plan<double>
99  {
100  typedef double scalar_type;
101  typedef fftw_complex complex_type;
102  std::shared_ptr<fftw_plan_s> m_plan;
103  fftw_plan() = default;
104 
105  void set_plan(::fftw_plan p) { m_plan.reset(p, fftw_destroy_plan); }
106  inline
107  void fwd(complex_type * dst,complex_type * src,int nfft) {
108  if (m_plan==NULL) set_plan(fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
109  fftw_execute_dft( m_plan.get(), src,dst);
110  }
111  inline
112  void inv(complex_type * dst,complex_type * src,int nfft) {
113  if (m_plan==NULL) set_plan(fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
114  fftw_execute_dft( m_plan.get(), src,dst);
115  }
116  inline
117  void fwd(complex_type * dst,scalar_type * src,int nfft) {
118  if (m_plan==NULL) set_plan(fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
119  fftw_execute_dft_r2c( m_plan.get(),src,dst);
120  }
121  inline
122  void inv(scalar_type * dst,complex_type * src,int nfft) {
123  if (m_plan==NULL)
124  set_plan(fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
125  fftw_execute_dft_c2r( m_plan.get(), src,dst);
126  }
127  inline
128  void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
129  if (m_plan==NULL) set_plan(fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
130  fftw_execute_dft( m_plan.get(), src,dst);
131  }
132  inline
133  void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
134  if (m_plan==NULL) set_plan(fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
135  fftw_execute_dft( m_plan.get(), src,dst);
136  }
137  };
138  template <>
139  struct fftw_plan<long double>
140  {
141  typedef long double scalar_type;
142  typedef fftwl_complex complex_type;
143  std::shared_ptr<fftwl_plan_s> m_plan;
144  fftw_plan() = default;
145 
146  void set_plan(fftwl_plan p) { m_plan.reset(p, fftwl_destroy_plan); }
147  inline
148  void fwd(complex_type * dst,complex_type * src,int nfft) {
149  if (m_plan==NULL) set_plan(fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
150  fftwl_execute_dft( m_plan.get(), src,dst);
151  }
152  inline
153  void inv(complex_type * dst,complex_type * src,int nfft) {
154  if (m_plan==NULL) set_plan(fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
155  fftwl_execute_dft( m_plan.get(), src,dst);
156  }
157  inline
158  void fwd(complex_type * dst,scalar_type * src,int nfft) {
159  if (m_plan==NULL) set_plan(fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
160  fftwl_execute_dft_r2c( m_plan.get(),src,dst);
161  }
162  inline
163  void inv(scalar_type * dst,complex_type * src,int nfft) {
164  if (m_plan==NULL)
165  set_plan(fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
166  fftwl_execute_dft_c2r( m_plan.get(), src,dst);
167  }
168  inline
169  void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
170  if (m_plan==NULL) set_plan(fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
171  fftwl_execute_dft( m_plan.get(), src,dst);
172  }
173  inline
174  void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
175  if (m_plan==NULL) set_plan(fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
176  fftwl_execute_dft( m_plan.get(), src,dst);
177  }
178  };
179 
180  template <typename Scalar_>
181  struct fftw_impl
182  {
183  typedef Scalar_ Scalar;
184  typedef std::complex<Scalar> Complex;
185 
186  inline
187  void clear()
188  {
189  m_plans.clear();
190  }
191 
192  // complex-to-complex forward FFT
193  inline
194  void fwd( Complex * dst,const Complex *src,int nfft)
195  {
196  get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src),nfft );
197  }
198 
199  // real-to-complex forward FFT
200  inline
201  void fwd( Complex * dst,const Scalar * src,int nfft)
202  {
203  get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src) ,nfft);
204  }
205 
206  // 2-d complex-to-complex
207  inline
208  void fwd2(Complex * dst, const Complex * src, int n0,int n1)
209  {
210  get_plan(n0,n1,false,dst,src).fwd2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
211  }
212 
213  // inverse complex-to-complex
214  inline
215  void inv(Complex * dst,const Complex *src,int nfft)
216  {
217  get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
218  }
219 
220  // half-complex to scalar
221  inline
222  void inv( Scalar * dst,const Complex * src,int nfft)
223  {
224  get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
225  }
226 
227  // 2-d complex-to-complex
228  inline
229  void inv2(Complex * dst, const Complex * src, int n0,int n1)
230  {
231  get_plan(n0,n1,true,dst,src).inv2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
232  }
233 
234 
235  protected:
236  typedef fftw_plan<Scalar> PlanData;
237 
238  typedef Eigen::numext::int64_t int64_t;
239 
240  typedef std::map<int64_t,PlanData> PlanMap;
241 
242  PlanMap m_plans;
243 
244  inline
245  PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
246  {
247  bool inplace = (dst==src);
248  bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
249  int64_t key = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1;
250  return m_plans[key];
251  }
252 
253  inline
254  PlanData & get_plan(int n0,int n1,bool inverse,void * dst,const void * src)
255  {
256  bool inplace = (dst==src);
257  bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
258  int64_t key = ( ( (((int64_t)n0) << 30)|(n1<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1;
259  return m_plans[key];
260  }
261  };
262 
263 } // end namespace internal
264 
265 } // end namespace Eigen
p
float * p
Eigen::internal::fftw_cast
T * fftw_cast(const T *p)
Definition: ei_fftw_impl.h:28
Eigen::numext::int64_t
std::int64_t int64_t
Eigen
: TensorContractionSycl.h, provides various tensor contraction kernel for SYCL backend
Eigen::inverse
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_inverse_op< typename Derived::Scalar >, const Derived > inverse(const Eigen::ArrayBase< Derived > &x)
Generated on Tue Jul 18 2023 15:28:44 for Eigen-unsupported by Doxygen 1.9.1 and the doxyYoda theme. This is an unofficial themed variant maintained by Rohit Goswami, based on the official documentation.