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 namespace Eigen {
11 
12 namespace 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;
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:
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: */
void inv(complex_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:65
void inv(scalar_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:118
void inv(scalar_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:75
void fwd(Complex *dst, const Complex *src, int nfft)
Definition: ei_fftw_impl.h:190
USING_NAMESPACE_ACADO typedef TaylorVariable< Interval > T
void fwd2(complex_type *dst, complex_type *src, int n0, int n1)
Definition: ei_fftw_impl.h:124
void inv2(complex_type *dst, complex_type *src, int n0, int n1)
Definition: ei_fftw_impl.h:87
void inv(Complex *dst, const Complex *src, int nfft)
Definition: ei_fftw_impl.h:211
fftw_plan< Scalar > PlanData
Definition: ei_fftw_impl.h:232
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
void fwd(Complex *dst, const Scalar *src, int nfft)
Definition: ei_fftw_impl.h:197
void inv2(complex_type *dst, complex_type *src, int n0, int n1)
Definition: ei_fftw_impl.h:129
void fwd(complex_type *dst, scalar_type *src, int nfft)
Definition: ei_fftw_impl.h:113
std::complex< Scalar > Complex
Definition: ei_fftw_impl.h:180
const CwiseUnaryOp< internal::scalar_inverse_op< Scalar >, const Derived > inverse() const
void inv(scalar_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:159
void inv(Scalar *dst, const Complex *src, int nfft)
Definition: ei_fftw_impl.h:218
void inv(complex_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:149
void inv2(Complex *dst, const Complex *src, int n0, int n1)
Definition: ei_fftw_impl.h:225
void fwd2(complex_type *dst, complex_type *src, int n0, int n1)
Definition: ei_fftw_impl.h:82
void fwd(complex_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:60
void fwd2(Complex *dst, const Complex *src, int n0, int n1)
Definition: ei_fftw_impl.h:204
void inv2(complex_type *dst, complex_type *src, int n0, int n1)
Definition: ei_fftw_impl.h:170
void fwd(complex_type *dst, scalar_type *src, int nfft)
Definition: ei_fftw_impl.h:70
std::map< int64_t, PlanData > PlanMap
Definition: ei_fftw_impl.h:234
PlanData & get_plan(int n0, int n1, bool inverse, void *dst, const void *src)
Definition: ei_fftw_impl.h:248
void fwd2(complex_type *dst, complex_type *src, int n0, int n1)
Definition: ei_fftw_impl.h:165
T * fftw_cast(const T *p)
Definition: ei_fftw_impl.h:24
void fwd(complex_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:103
void inv(complex_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:108
void fwd(complex_type *dst, scalar_type *src, int nfft)
Definition: ei_fftw_impl.h:154
PlanData & get_plan(int nfft, bool inverse, void *dst, const void *src)
Definition: ei_fftw_impl.h:239
void fwd(complex_type *dst, complex_type *src, int nfft)
Definition: ei_fftw_impl.h:144


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:32