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