From 00b75375e73b229466dcd54825aad52d5bb4e9ea Mon Sep 17 00:00:00 2001 From: Guoqiang QI Date: Wed, 11 May 2022 17:44:22 +0000 Subject: [PATCH] Adding PocketFFT support in FFT module since kissfft has some flaw in accuracy and performance --- test/main.h | 2 +- unsupported/Eigen/FFT | 32 +- unsupported/Eigen/src/FFT/ei_pocketfft_impl.h | 69 +++++ unsupported/test/CMakeLists.txt | 11 + unsupported/test/FFT.cpp | 4 +- unsupported/test/FFTW.cpp | 264 +---------------- unsupported/test/fft_test_shared.h | 279 ++++++++++++++++++ unsupported/test/pocketfft.cpp | 2 + 8 files changed, 390 insertions(+), 273 deletions(-) create mode 100644 unsupported/Eigen/src/FFT/ei_pocketfft_impl.h create mode 100644 unsupported/test/fft_test_shared.h create mode 100644 unsupported/test/pocketfft.cpp diff --git a/test/main.h b/test/main.h index 8eead7cd7..884b6c34d 100644 --- a/test/main.h +++ b/test/main.h @@ -87,7 +87,7 @@ // protected by parenthesis against macro expansion, the min()/max() macros // are defined here and any not-parenthesized min/max call will cause a // compiler error. -#if !defined(__HIPCC__) && !defined(EIGEN_USE_SYCL) +#if !defined(__HIPCC__) && !defined(EIGEN_USE_SYCL) && !defined(EIGEN_POCKETFFT_DEFAULT) // // HIP header files include the following files // diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index c9d693897..59bb7165e 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -29,10 +29,19 @@ * The default implementation is based on kissfft. It is a small, free, and * reasonably efficient default. * - * There are currently two implementation backend: + * There are currently four implementation backend: * + * - kissfft(https://github.com/mborgerding/kissfft) : Simple and not so fast, BSD-3-Clause. + * It is a mixed-radix Fast Fourier Transform based up on the principle, "Keep It Simple, Stupid." + * Notice that:kissfft fails to handle "atypically-sized" inputs(i.e., sizes with large factors),a workaround is using fftw or pocketfft. * - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size. * - MKL (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form. + * - pocketfft (https://gitlab.mpcdf.mpg.de/mtr/pocketfft) : faster than kissfft, BSD 3-clause. + * It is a heavily modified implementation of FFTPack, with the following advantages: + * 1.strictly C++11 compliant + * 2.more accurate twiddle factor computation + * 3.very fast plan generation + * 4.worst case complexity for transform sizes with large prime factors is N*log(N), because Bluestein's algorithm is used for these cases. * * \section FFTDesign Design * @@ -85,9 +94,16 @@ namespace Eigen { template struct default_fft_impl : public internal::imklfft_impl {}; } -#else +#elif defined EIGEN_POCKETFFT_DEFAULT +// internal::pocketfft_impl: a heavily modified implementation of FFTPack, with many advantages. +# include +# include"src/FFT/ei_pocketfft_impl.h" + namespace Eigen { + template + struct default_fft_impl : public internal::pocketfft_impl {}; + } +#else // internal::kissfft_impl: small, free, reasonably efficient default, derived from kissfft -// # include "src/FFT/ei_kissfft_impl.h" namespace Eigen { template @@ -195,13 +211,13 @@ class FFT m_impl.fwd(dst,src,static_cast(nfft)); } - /* +#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT inline void fwd2(Complex * dst, const Complex * src, int n0,int n1) { m_impl.fwd2(dst,src,n0,n1); } - */ +#endif template inline @@ -354,8 +370,7 @@ class FFT } - /* - // TODO: multi-dimensional FFTs +#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT inline void inv2(Complex * dst, const Complex * src, int n0,int n1) { @@ -363,7 +378,8 @@ class FFT if ( HasFlag( Unscaled ) == false) scale(dst,1./(n0*n1),n0*n1); } - */ +#endif + inline impl_type & impl() {return m_impl;} diff --git a/unsupported/Eigen/src/FFT/ei_pocketfft_impl.h b/unsupported/Eigen/src/FFT/ei_pocketfft_impl.h new file mode 100644 index 000000000..f2da890c7 --- /dev/null +++ b/unsupported/Eigen/src/FFT/ei_pocketfft_impl.h @@ -0,0 +1,69 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +using namespace pocketfft; +using namespace pocketfft::detail; + +namespace Eigen { + +namespace internal { + +template +struct pocketfft_impl +{ + typedef _Scalar Scalar; + typedef std::complex Complex; + + inline void clear() {} + + inline void fwd(Complex* dst, const Scalar* src, int nfft){ + const shape_t shape_{ static_cast(nfft) }; + const shape_t axes_{ 0 }; + const stride_t stride_in{ sizeof(Scalar) }; + const stride_t stride_out{ sizeof(Complex) }; + r2c(shape_, stride_in, stride_out, axes_, FORWARD, src, dst, static_cast(1)); + } + + inline void fwd(Complex* dst, const Complex* src, int nfft){ + const shape_t shape_{ static_cast(nfft) }; + const shape_t axes_{ 0 }; + const stride_t stride_{ sizeof(Complex) }; + c2c(shape_, stride_, stride_, axes_, FORWARD, src, dst, static_cast(1)); + } + + inline void inv(Scalar* dst, const Complex* src, int nfft){ + const shape_t shape_{ static_cast(nfft) }; + const shape_t axes_{ 0 }; + const stride_t stride_in{ sizeof(Complex) }; + const stride_t stride_out{ sizeof(Scalar) }; + c2r(shape_, stride_in, stride_out, axes_, BACKWARD, src, dst, static_cast(1)); + } + + inline void inv(Complex* dst, const Complex* src, int nfft){ + const shape_t shape_{ static_cast(nfft) }; + const shape_t axes_{ 0 }; + const stride_t stride_{ sizeof(Complex) }; + c2c(shape_, stride_, stride_, axes_, BACKWARD, src, dst, static_cast(1)); + } + + inline void fwd2(Complex* dst, const Complex* src, int nfft0, int nfft1){ + const shape_t shape_{ static_cast(nfft0), static_cast(nfft1) }; + const shape_t axes_{ 0, 1 }; + const stride_t stride_{ static_cast(sizeof(Complex)*nfft1), static_cast(sizeof(Complex)) }; + c2c(shape_, stride_, stride_, axes_, FORWARD, src, dst, static_cast(1)); + } + + inline void inv2(Complex* dst, const Complex* src, int nfft0, int nfft1){ + const shape_t shape_{ static_cast(nfft0), static_cast(nfft1) }; + const shape_t axes_{ 0, 1 }; + const stride_t stride_{ static_cast(sizeof(Complex)*nfft1), static_cast(sizeof(Complex)) }; + c2c(shape_, stride_, stride_, axes_, BACKWARD, src, dst, static_cast(1)); + } +}; + +} // namespace internal +} // namespace Eigen diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index b7661e70c..21d8c5eb5 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -77,6 +77,17 @@ else() ei_add_property(EIGEN_MISSING_BACKENDS "fftw, ") endif() +find_path(POCKETFFT pocketfft_hdronly.h) +if(POCKETFFT) + if(EIGEN_TEST_CXX11) + ei_add_property(EIGEN_TESTED_BACKENDS "pocketfft, ") + include_directories( ${POCKETFFT} ) + ei_add_test(pocketfft "-pthread" "${CMAKE_THREAD_LIBS_INIT}" "-DEIGEN_POCKETFFT_DEFAULT" ) + endif() +else() + ei_add_property(EIGEN_MISSING_BACKENDS "pocketfft, ") +endif() + option(EIGEN_TEST_OPENGL "Enable OpenGL support in unit tests" OFF) if(EIGEN_TEST_OPENGL) find_package(OpenGL) diff --git a/unsupported/test/FFT.cpp b/unsupported/test/FFT.cpp index 45c87f5a7..f85461cc5 100644 --- a/unsupported/test/FFT.cpp +++ b/unsupported/test/FFT.cpp @@ -1,2 +1,2 @@ -#define test_FFTW test_FFT -#include "FFTW.cpp" +#define EIGEN_FFT_DEFAULT 1 +#include "fft_test_shared.h" diff --git a/unsupported/test/FFTW.cpp b/unsupported/test/FFTW.cpp index cfe559ebd..d69867c7e 100644 --- a/unsupported/test/FFTW.cpp +++ b/unsupported/test/FFTW.cpp @@ -1,262 +1,2 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Mark Borgerding mark a borgerding net -// -// This Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "main.h" -#include - -template -std::complex RandomCpx() { return std::complex( (T)(rand()/(T)RAND_MAX - .5), (T)(rand()/(T)RAND_MAX - .5) ); } - -using namespace std; -using namespace Eigen; - - -template < typename T> -complex promote(complex x) { return complex((long double)x.real(),(long double)x.imag()); } - -complex promote(float x) { return complex((long double)x); } -complex promote(double x) { return complex((long double)x); } -complex promote(long double x) { return complex((long double)x); } - - - template - long double fft_rmse( const VT1 & fftbuf,const VT2 & timebuf) - { - long double totalpower=0; - long double difpower=0; - long double pi = acos((long double)-1 ); - for (size_t k0=0;k0<(size_t)fftbuf.size();++k0) { - complex acc = 0; - long double phinc = (long double)(-2.)*k0* pi / timebuf.size(); - for (size_t k1=0;k1<(size_t)timebuf.size();++k1) { - acc += promote( timebuf[k1] ) * exp( complex(0,k1*phinc) ); - } - totalpower += numext::abs2(acc); - complex x = promote(fftbuf[k0]); - complex dif = acc - x; - difpower += numext::abs2(dif); - //cerr << k0 << "\t" << acc << "\t" << x << "\t" << sqrt(numext::abs2(dif)) << endl; - } - cerr << "rmse:" << sqrt(difpower/totalpower) << endl; - return sqrt(difpower/totalpower); - } - - template - long double dif_rmse( const VT1 buf1,const VT2 buf2) - { - long double totalpower=0; - long double difpower=0; - size_t n = (min)( buf1.size(),buf2.size() ); - for (size_t k=0;k struct VectorType; - -template struct VectorType -{ - typedef vector type; -}; - -template struct VectorType -{ - typedef Matrix type; -}; - -template -void test_scalar_generic(int nfft) -{ - typedef typename FFT::Complex Complex; - typedef typename FFT::Scalar Scalar; - typedef typename VectorType::type ScalarVector; - typedef typename VectorType::type ComplexVector; - - FFT fft; - ScalarVector tbuf(nfft); - ComplexVector freqBuf; - for (int k=0;k>1)+1) ); - VERIFY( T(fft_rmse(freqBuf,tbuf)) < test_precision() );// gross check - - fft.ClearFlag(fft.HalfSpectrum ); - fft.fwd( freqBuf,tbuf); - VERIFY( (size_t)freqBuf.size() == (size_t)nfft); - VERIFY( T(fft_rmse(freqBuf,tbuf)) < test_precision() );// gross check - - if (nfft&1) - return; // odd FFTs get the wrong size inverse FFT - - ScalarVector tbuf2; - fft.inv( tbuf2 , freqBuf); - VERIFY( T(dif_rmse(tbuf,tbuf2)) < test_precision() );// gross check - - - // verify that the Unscaled flag takes effect - ScalarVector tbuf3; - fft.SetFlag(fft.Unscaled); - - fft.inv( tbuf3 , freqBuf); - - for (int k=0;k " << (tbuf3[i] - tbuf[i] ) << endl; - - VERIFY( T(dif_rmse(tbuf,tbuf3)) < test_precision() );// gross check - - // verify that ClearFlag works - fft.ClearFlag(fft.Unscaled); - fft.inv( tbuf2 , freqBuf); - VERIFY( T(dif_rmse(tbuf,tbuf2)) < test_precision() );// gross check -} - -template -void test_scalar(int nfft) -{ - test_scalar_generic(nfft); - //test_scalar_generic(nfft); -} - - -template -void test_complex_generic(int nfft) -{ - typedef typename FFT::Complex Complex; - typedef typename VectorType::type ComplexVector; - - FFT fft; - - ComplexVector inbuf(nfft); - ComplexVector outbuf; - ComplexVector buf3; - for (int k=0;k() );// gross check - fft.inv( buf3 , outbuf); - - VERIFY( T(dif_rmse(inbuf,buf3)) < test_precision() );// gross check - - // verify that the Unscaled flag takes effect - ComplexVector buf4; - fft.SetFlag(fft.Unscaled); - fft.inv( buf4 , outbuf); - for (int k=0;k() );// gross check - - // verify that ClearFlag works - fft.ClearFlag(fft.Unscaled); - fft.inv( buf3 , outbuf); - VERIFY( T(dif_rmse(inbuf,buf3)) < test_precision() );// gross check -} - -template -void test_complex(int nfft) -{ - test_complex_generic(nfft); - test_complex_generic(nfft); -} -/* -template -void test_complex2d() -{ - typedef typename Eigen::FFT::Complex Complex; - FFT fft; - Eigen::Matrix src,src2,dst,dst2; - - src = Eigen::Matrix::Random(); - //src = Eigen::Matrix::Identity(); - - for (int k=0;k tmpOut; - fft.fwd( tmpOut,src.col(k) ); - dst2.col(k) = tmpOut; - } - - for (int k=0;k tmpOut; - fft.fwd( tmpOut, dst2.row(k) ); - dst2.row(k) = tmpOut; - } - - fft.fwd2(dst.data(),src.data(),ncols,nrows); - fft.inv2(src2.data(),dst.data(),ncols,nrows); - VERIFY( (src-src2).norm() < test_precision() ); - VERIFY( (dst-dst2).norm() < test_precision() ); -} -*/ - - -void test_return_by_value(int len) -{ - VectorXf in; - VectorXf in1; - in.setRandom( len ); - VectorXcf out1,out2; - FFT fft; - - fft.SetFlag(fft.HalfSpectrum ); - - fft.fwd(out1,in); - out2 = fft.fwd(in); - VERIFY( (out1-out2).norm() < test_precision() ); - in1 = fft.inv(out1); - VERIFY( (in1-in).norm() < test_precision() ); -} - -EIGEN_DECLARE_TEST(FFTW) -{ - CALL_SUBTEST( test_return_by_value(32) ); - //CALL_SUBTEST( ( test_complex2d () ) ); CALL_SUBTEST( ( test_complex2d () ) ); - //CALL_SUBTEST( ( test_complex2d () ) ); - CALL_SUBTEST( test_complex(32) ); CALL_SUBTEST( test_complex(32) ); - CALL_SUBTEST( test_complex(256) ); CALL_SUBTEST( test_complex(256) ); - CALL_SUBTEST( test_complex(3*8) ); CALL_SUBTEST( test_complex(3*8) ); - CALL_SUBTEST( test_complex(5*32) ); CALL_SUBTEST( test_complex(5*32) ); - CALL_SUBTEST( test_complex(2*3*4) ); CALL_SUBTEST( test_complex(2*3*4) ); - CALL_SUBTEST( test_complex(2*3*4*5) ); CALL_SUBTEST( test_complex(2*3*4*5) ); - CALL_SUBTEST( test_complex(2*3*4*5*7) ); CALL_SUBTEST( test_complex(2*3*4*5*7) ); - - CALL_SUBTEST( test_scalar(32) ); CALL_SUBTEST( test_scalar(32) ); - CALL_SUBTEST( test_scalar(45) ); CALL_SUBTEST( test_scalar(45) ); - CALL_SUBTEST( test_scalar(50) ); CALL_SUBTEST( test_scalar(50) ); - CALL_SUBTEST( test_scalar(256) ); CALL_SUBTEST( test_scalar(256) ); - CALL_SUBTEST( test_scalar(2*3*4*5*7) ); CALL_SUBTEST( test_scalar(2*3*4*5*7) ); - - #ifdef EIGEN_HAS_FFTWL - CALL_SUBTEST( test_complex(32) ); - CALL_SUBTEST( test_complex(256) ); - CALL_SUBTEST( test_complex(3*8) ); - CALL_SUBTEST( test_complex(5*32) ); - CALL_SUBTEST( test_complex(2*3*4) ); - CALL_SUBTEST( test_complex(2*3*4*5) ); - CALL_SUBTEST( test_complex(2*3*4*5*7) ); - - CALL_SUBTEST( test_scalar(32) ); - CALL_SUBTEST( test_scalar(45) ); - CALL_SUBTEST( test_scalar(50) ); - CALL_SUBTEST( test_scalar(256) ); - CALL_SUBTEST( test_scalar(2*3*4*5*7) ); - #endif -} +#define EIGEN_FFTW_DEFAULT 1 +#include "fft_test_shared.h" diff --git a/unsupported/test/fft_test_shared.h b/unsupported/test/fft_test_shared.h new file mode 100644 index 000000000..06cf365f6 --- /dev/null +++ b/unsupported/test/fft_test_shared.h @@ -0,0 +1,279 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Mark Borgerding mark a borgerding net +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include + +template +std::complex RandomCpx() { return std::complex( (T)(rand()/(T)RAND_MAX - .5), (T)(rand()/(T)RAND_MAX - .5) ); } + +using namespace std; +using namespace Eigen; + + +template < typename T> +complex promote(complex x) { return complex((long double)x.real(),(long double)x.imag()); } + +complex promote(float x) { return complex((long double)x); } +complex promote(double x) { return complex((long double)x); } +complex promote(long double x) { return complex((long double)x); } + + + template + long double fft_rmse( const VT1 & fftbuf,const VT2 & timebuf) + { + long double totalpower=0; + long double difpower=0; + long double pi = acos((long double)-1 ); + for (size_t k0=0;k0<(size_t)fftbuf.size();++k0) { + complex acc = 0; + long double phinc = (long double)(-2.)*k0* pi / timebuf.size(); + for (size_t k1=0;k1<(size_t)timebuf.size();++k1) { + acc += promote( timebuf[k1] ) * exp( complex(0,k1*phinc) ); + } + totalpower += numext::abs2(acc); + complex x = promote(fftbuf[k0]); + complex dif = acc - x; + difpower += numext::abs2(dif); + //cerr << k0 << "\t" << acc << "\t" << x << "\t" << sqrt(numext::abs2(dif)) << endl; + } + // cerr << "rmse:" << sqrt(difpower/totalpower) << endl; + return sqrt(difpower/totalpower); + } + + template + long double dif_rmse( const VT1 buf1,const VT2 buf2) + { + long double totalpower=0; + long double difpower=0; + size_t n = (min)( buf1.size(),buf2.size() ); + for (size_t k=0;k struct VectorType; + +template struct VectorType +{ + typedef vector type; +}; + +template struct VectorType +{ + typedef Matrix type; +}; + +template +void test_scalar_generic(int nfft) +{ + typedef typename FFT::Complex Complex; + typedef typename FFT::Scalar Scalar; + typedef typename VectorType::type ScalarVector; + typedef typename VectorType::type ComplexVector; + + FFT fft; + ScalarVector tbuf(nfft); + ComplexVector freqBuf; + for (int k=0;k>1)+1) ); + VERIFY( T(fft_rmse(freqBuf,tbuf)) < test_precision() );// gross check + + fft.ClearFlag(fft.HalfSpectrum ); + fft.fwd( freqBuf,tbuf); + VERIFY( (size_t)freqBuf.size() == (size_t)nfft); + VERIFY( T(fft_rmse(freqBuf,tbuf)) < test_precision() );// gross check + + if (nfft&1) + return; // odd FFTs get the wrong size inverse FFT + + ScalarVector tbuf2; + fft.inv( tbuf2 , freqBuf); + VERIFY( T(dif_rmse(tbuf,tbuf2)) < test_precision() );// gross check + + + // verify that the Unscaled flag takes effect + ScalarVector tbuf3; + fft.SetFlag(fft.Unscaled); + + fft.inv( tbuf3 , freqBuf); + + for (int k=0;k " << (tbuf3[i] - tbuf[i] ) << endl; + + VERIFY( T(dif_rmse(tbuf,tbuf3)) < test_precision() );// gross check + + // verify that ClearFlag works + fft.ClearFlag(fft.Unscaled); + fft.inv( tbuf2 , freqBuf); + VERIFY( T(dif_rmse(tbuf,tbuf2)) < test_precision() );// gross check +} + +template +void test_scalar(int nfft) +{ + test_scalar_generic(nfft); + //test_scalar_generic(nfft); +} + + +template +void test_complex_generic(int nfft) +{ + typedef typename FFT::Complex Complex; + typedef typename VectorType::type ComplexVector; + + FFT fft; + + ComplexVector inbuf(nfft); + ComplexVector outbuf; + ComplexVector buf3; + for (int k=0;k() );// gross check + fft.inv( buf3 , outbuf); + + VERIFY( T(dif_rmse(inbuf,buf3)) < test_precision() );// gross check + + // verify that the Unscaled flag takes effect + ComplexVector buf4; + fft.SetFlag(fft.Unscaled); + fft.inv( buf4 , outbuf); + for (int k=0;k() );// gross check + + // verify that ClearFlag works + fft.ClearFlag(fft.Unscaled); + fft.inv( buf3 , outbuf); + VERIFY( T(dif_rmse(inbuf,buf3)) < test_precision() );// gross check +} + +template +void test_complex(int nfft) +{ + test_complex_generic(nfft); + test_complex_generic(nfft); +} + +template +void test_complex2d() +{ + typedef typename Eigen::FFT::Complex Complex; + FFT fft; + Eigen::Matrix src,src2,dst,dst2; + + src = Eigen::Matrix::Random(); + //src = Eigen::Matrix::Identity(); + + for (int k=0;k tmpOut; + fft.fwd( tmpOut,src.col(k) ); + dst2.col(k) = tmpOut; + } + + for (int k=0;k tmpOut; + fft.fwd( tmpOut, dst2.row(k) ); + dst2.row(k) = tmpOut; + } + + fft.fwd2(dst.data(),src.data(),ncols,nrows); + fft.inv2(src2.data(),dst.data(),ncols,nrows); + VERIFY( (src-src2).norm() < test_precision() ); + VERIFY( (dst-dst2).norm() < test_precision() ); +} + +void test_return_by_value(int len) +{ + VectorXf in; + VectorXf in1; + in.setRandom( len ); + VectorXcf out1,out2; + FFT fft; + + fft.SetFlag(fft.HalfSpectrum ); + + fft.fwd(out1,in); + out2 = fft.fwd(in); + VERIFY( (out1-out2).norm() < test_precision() ); + in1 = fft.inv(out1); + VERIFY( (in1-in).norm() < test_precision() ); +} + +EIGEN_DECLARE_TEST(FFTW) +{ + CALL_SUBTEST( test_return_by_value(32) ); + CALL_SUBTEST( test_complex(32) ); CALL_SUBTEST( test_complex(32) ); + CALL_SUBTEST( test_complex(256) ); CALL_SUBTEST( test_complex(256) ); + CALL_SUBTEST( test_complex(3*8) ); CALL_SUBTEST( test_complex(3*8) ); + CALL_SUBTEST( test_complex(5*32) ); CALL_SUBTEST( test_complex(5*32) ); + CALL_SUBTEST( test_complex(2*3*4) ); CALL_SUBTEST( test_complex(2*3*4) ); + CALL_SUBTEST( test_complex(2*3*4*5) ); CALL_SUBTEST( test_complex(2*3*4*5) ); + CALL_SUBTEST( test_complex(2*3*4*5*7) ); CALL_SUBTEST( test_complex(2*3*4*5*7) ); + + CALL_SUBTEST( test_scalar(32) ); CALL_SUBTEST( test_scalar(32) ); + CALL_SUBTEST( test_scalar(45) ); CALL_SUBTEST( test_scalar(45) ); + CALL_SUBTEST( test_scalar(50) ); CALL_SUBTEST( test_scalar(50) ); + CALL_SUBTEST( test_scalar(256) ); CALL_SUBTEST( test_scalar(256) ); + CALL_SUBTEST( test_scalar(2*3*4*5*7) ); CALL_SUBTEST( test_scalar(2*3*4*5*7) ); + + #if defined EIGEN_HAS_FFTWL || defined EIGEN_POCKETFFT_DEFAULT + CALL_SUBTEST( test_complex(32) ); + CALL_SUBTEST( test_complex(256) ); + CALL_SUBTEST( test_complex(3*8) ); + CALL_SUBTEST( test_complex(5*32) ); + CALL_SUBTEST( test_complex(2*3*4) ); + CALL_SUBTEST( test_complex(2*3*4*5) ); + CALL_SUBTEST( test_complex(2*3*4*5*7) ); + + CALL_SUBTEST( test_scalar(32) ); + CALL_SUBTEST( test_scalar(45) ); + CALL_SUBTEST( test_scalar(50) ); + CALL_SUBTEST( test_scalar(256) ); + CALL_SUBTEST( test_scalar(2*3*4*5*7) ); + + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + // fail to build since Eigen limit the stack allocation size,too big here. + // CALL_SUBTEST( ( test_complex2d () ) ); + + #endif + + #if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + #endif + +} diff --git a/unsupported/test/pocketfft.cpp b/unsupported/test/pocketfft.cpp new file mode 100644 index 000000000..5e8a0b6ae --- /dev/null +++ b/unsupported/test/pocketfft.cpp @@ -0,0 +1,2 @@ +#define EIGEN_POCKETFFT_DEFAULT 1 +#include "fft_test_shared.h"