ei_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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 namespace Eigen {
26 
27 namespace internal {
28 
29  // FFTW uses non-const arguments
30  // so we must use ugly const_cast calls for all the args it uses
31  //
32  // This should be safe as long as
33  // 1. we use FFTW_ESTIMATE for all our planning
34  // see the FFTW docs section 4.3.2 "Planner Flags"
35  // 2. fftw_complex is compatible with std::complex
36  // This assumes std::complex<T> layout is array of size 2 with real,imag
37  template <typename T>
38  inline
39  T * fftw_cast(const T* p)
40  {
41  return const_cast<T*>( p);
42  }
43 
44  inline
45  fftw_complex * fftw_cast( const std::complex<double> * p)
46  {
47  return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) );
48  }
49 
50  inline
51  fftwf_complex * fftw_cast( const std::complex<float> * p)
52  {
53  return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) );
54  }
55 
56  inline
57  fftwl_complex * fftw_cast( const std::complex<long double> * p)
58  {
59  return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) );
60  }
61 
62  template <typename T>
63  struct fftw_plan {};
64 
65  template <>
66  struct fftw_plan<float>
67  {
68  typedef float scalar_type;
69  typedef fftwf_complex complex_type;
70  fftwf_plan m_plan;
71  fftw_plan() :m_plan(NULL) {}
72  ~fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
73 
74  inline
75  void fwd(complex_type * dst,complex_type * src,int nfft) {
76  if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
77  fftwf_execute_dft( m_plan, src,dst);
78  }
79  inline
80  void inv(complex_type * dst,complex_type * src,int nfft) {
81  if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
82  fftwf_execute_dft( m_plan, src,dst);
83  }
84  inline
85  void fwd(complex_type * dst,scalar_type * src,int nfft) {
86  if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
87  fftwf_execute_dft_r2c( m_plan,src,dst);
88  }
89  inline
90  void inv(scalar_type * dst,complex_type * src,int nfft) {
91  if (m_plan==NULL)
92  m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
93  fftwf_execute_dft_c2r( m_plan, src,dst);
94  }
95 
96  inline
97  void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
98  if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
99  fftwf_execute_dft( m_plan, src,dst);
100  }
101  inline
102  void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
103  if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
104  fftwf_execute_dft( m_plan, src,dst);
105  }
106 
107  };
108  template <>
109  struct fftw_plan<double>
110  {
111  typedef double scalar_type;
112  typedef fftw_complex complex_type;
113  ::fftw_plan m_plan;
114  fftw_plan() :m_plan(NULL) {}
115  ~fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
116 
117  inline
118  void fwd(complex_type * dst,complex_type * src,int nfft) {
119  if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
120  fftw_execute_dft( m_plan, src,dst);
121  }
122  inline
123  void inv(complex_type * dst,complex_type * src,int nfft) {
124  if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
125  fftw_execute_dft( m_plan, src,dst);
126  }
127  inline
128  void fwd(complex_type * dst,scalar_type * src,int nfft) {
129  if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
130  fftw_execute_dft_r2c( m_plan,src,dst);
131  }
132  inline
133  void inv(scalar_type * dst,complex_type * src,int nfft) {
134  if (m_plan==NULL)
135  m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
136  fftw_execute_dft_c2r( m_plan, src,dst);
137  }
138  inline
139  void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
140  if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
141  fftw_execute_dft( m_plan, src,dst);
142  }
143  inline
144  void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
145  if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
146  fftw_execute_dft( m_plan, src,dst);
147  }
148  };
149  template <>
150  struct fftw_plan<long double>
151  {
152  typedef long double scalar_type;
153  typedef fftwl_complex complex_type;
154  fftwl_plan m_plan;
155  fftw_plan() :m_plan(NULL) {}
156  ~fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
157 
158  inline
159  void fwd(complex_type * dst,complex_type * src,int nfft) {
160  if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
161  fftwl_execute_dft( m_plan, src,dst);
162  }
163  inline
164  void inv(complex_type * dst,complex_type * src,int nfft) {
165  if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
166  fftwl_execute_dft( m_plan, src,dst);
167  }
168  inline
169  void fwd(complex_type * dst,scalar_type * src,int nfft) {
170  if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
171  fftwl_execute_dft_r2c( m_plan,src,dst);
172  }
173  inline
174  void inv(scalar_type * dst,complex_type * src,int nfft) {
175  if (m_plan==NULL)
176  m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
177  fftwl_execute_dft_c2r( m_plan, src,dst);
178  }
179  inline
180  void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
181  if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
182  fftwl_execute_dft( m_plan, src,dst);
183  }
184  inline
185  void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
186  if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
187  fftwl_execute_dft( m_plan, src,dst);
188  }
189  };
190 
191  template <typename _Scalar>
192  struct fftw_impl
193  {
194  typedef _Scalar Scalar;
195  typedef std::complex<Scalar> Complex;
196 
197  inline
198  void clear()
199  {
200  m_plans.clear();
201  }
202 
203  // complex-to-complex forward FFT
204  inline
205  void fwd( Complex * dst,const Complex *src,int nfft)
206  {
207  get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src),nfft );
208  }
209 
210  // real-to-complex forward FFT
211  inline
212  void fwd( Complex * dst,const Scalar * src,int nfft)
213  {
214  get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src) ,nfft);
215  }
216 
217  // 2-d complex-to-complex
218  inline
219  void fwd2(Complex * dst, const Complex * src, int n0,int n1)
220  {
221  get_plan(n0,n1,false,dst,src).fwd2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
222  }
223 
224  // inverse complex-to-complex
225  inline
226  void inv(Complex * dst,const Complex *src,int nfft)
227  {
228  get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
229  }
230 
231  // half-complex to scalar
232  inline
233  void inv( Scalar * dst,const Complex * src,int nfft)
234  {
235  get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
236  }
237 
238  // 2-d complex-to-complex
239  inline
240  void inv2(Complex * dst, const Complex * src, int n0,int n1)
241  {
242  get_plan(n0,n1,true,dst,src).inv2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
243  }
244 
245 
246  protected:
247  typedef fftw_plan<Scalar> PlanData;
248 
249  typedef std::map<int64_t,PlanData> PlanMap;
250 
251  PlanMap m_plans;
252 
253  inline
254  PlanData & get_plan(int nfft,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 = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1;
259  return m_plans[key];
260  }
261 
262  inline
263  PlanData & get_plan(int n0,int n1,bool inverse,void * dst,const void * src)
264  {
265  bool inplace = (dst==src);
266  bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
267  int64_t key = ( ( (((int64_t)n0) << 30)|(n1<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1;
268  return m_plans[key];
269  }
270  };
271 
272 } // end namespace internal
273 
274 } // end namespace Eigen
275 
276 /* vim: set filetype=cpp et sw=2 ts=2 ai: */