diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 387c11385..0e8b49f75 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -190,6 +190,25 @@ template class MatrixBase * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */ int innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); } + /** Only plain matrices, not expressions may be resized; therefore the only useful resize method is + * Matrix::resize(). The present method only asserts that the new size equals the old size, and does + * nothing else. + */ + void resize(int size) + { + ei_assert(size == this->size() + && "MatrixBase::resize() does not actually allow to resize."); + } + /** Only plain matrices, not expressions may be resized; therefore the only useful resize method is + * Matrix::resize(). The present method only asserts that the new size equals the old size, and does + * nothing else. + */ + void resize(int rows, int cols) + { + ei_assert(rows == this->rows() && cols == this->cols() + && "MatrixBase::resize() does not actually allow to resize."); + } + #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal the plain matrix type corresponding to this expression. Note that is not necessarily * exactly the return type of eval(): in the case of plain matrices, the return type of eval() is a const diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h index 883f2d95e..6210bc91c 100644 --- a/Eigen/src/Core/util/StaticAssert.h +++ b/Eigen/src/Core/util/StaticAssert.h @@ -78,7 +78,8 @@ INVALID_MATRIX_TEMPLATE_PARAMETERS, BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER, THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX, - THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE + THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE, + THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES }; }; diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index dc7e85908..fafdb829b 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -36,7 +36,7 @@ #define DEFAULT_FFT_IMPL ei_fftw_impl #endif -// intel Math Kernel Library: fastest, commerical -- incompatible with Eigen in GPL form +// intel Math Kernel Library: fastest, commercial -- incompatible with Eigen in GPL form #ifdef _MKL_DFTI_H_ // mkl_dfti.h has been included, we can use MKL FFT routines // TODO // #include "src/FFT/ei_imkl_impl.h" @@ -70,6 +70,20 @@ class FFT fwd( &dst[0],&src[0],src.size() ); } + template + void fwd( MatrixBase & dst, const MatrixBase & src) + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) + EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,InputDerived) // size at compile-time + EIGEN_STATIC_ASSERT((ei_is_same_type::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit, + THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES) + dst.derived().resize( src.size() ); + fwd( &dst[0],&src[0],src.size() ); + } + template void inv( _Output * dst, const Complex * src, int nfft) { @@ -83,8 +97,24 @@ class FFT inv( &dst[0],&src[0],src.size() ); } + template + void inv( MatrixBase & dst, const MatrixBase & src) + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) + EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,OutputDerived) // size at compile-time + EIGEN_STATIC_ASSERT((ei_is_same_type::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + EIGEN_STATIC_ASSERT(int(OutputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit, + THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES) + dst.derived().resize( src.size() ); + inv( &dst[0],&src[0],src.size() ); + } + // TODO: multi-dimensional FFTs + // TODO: handle Eigen MatrixBase + // ---> i added fwd and inv specializations above + unit test, is this enough? (bjacob) traits_type & traits() {return m_traits;} private: diff --git a/unsupported/test/FFT.cpp b/unsupported/test/FFT.cpp index f0b9b68bf..cc68f3718 100644 --- a/unsupported/test/FFT.cpp +++ b/unsupported/test/FFT.cpp @@ -39,16 +39,16 @@ complex promote(double x) { return complex( x); } complex promote(long double x) { return complex( x); } - template - long double fft_rmse( const vector & fftbuf,const vector & timebuf) + template + long double fft_rmse( const VectorType1 & fftbuf,const VectorType2 & timebuf) { long double totalpower=0; long double difpower=0; cerr <<"idx\ttruth\t\tvalue\t|dif|=\n"; - for (size_t k0=0;k0 acc = 0; long double phinc = -2.*k0* M_PIl / timebuf.size(); - for (size_t k1=0;k1(0,k1*phinc) ); } totalpower += norm(acc); @@ -61,8 +61,8 @@ complex promote(long double x) { return complex( x); return sqrt(difpower/totalpower); } - template - long double dif_rmse( const vector buf1,const vector buf2) + template + long double dif_rmse( const VectorType1& buf1,const VectorType2& buf2) { long double totalpower=0; long double difpower=0; @@ -74,35 +74,59 @@ complex promote(long double x) { return complex( x); return sqrt(difpower/totalpower); } -template -void test_scalar(int nfft) +enum { StdVectorContainer, EigenVectorContainer }; + +template struct VectorType; + +template struct VectorType { - typedef typename Eigen::FFT::Complex Complex; - typedef typename Eigen::FFT::Scalar Scalar; + 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; - vector inbuf(nfft); - vector outbuf; + ScalarVector inbuf(nfft); + ComplexVector outbuf; for (int k=0;k() );// gross check - vector buf3; + ScalarVector buf3; fft.inv( buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision() );// gross check } -template -void test_complex(int nfft) +template +void test_scalar(int nfft) { - typedef typename Eigen::FFT::Complex Complex; + 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; - vector inbuf(nfft); - vector outbuf; - vector buf3; + ComplexVector inbuf(nfft); + ComplexVector outbuf; + ComplexVector buf3; for (int k=0;k() );// gross check } +template +void test_complex(int nfft) +{ + test_complex_generic(nfft); + test_complex_generic(nfft); +} + void test_FFT() { - CALL_SUBTEST( test_complex(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(256) ); - CALL_SUBTEST( test_complex(3*8) ); 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(5*32) ); - CALL_SUBTEST( test_complex(2*3*4) ); 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) ); - CALL_SUBTEST( test_complex(2*3*4*5*7) ); CALL_SUBTEST( test_complex(2*3*4*5*7) ); CALL_SUBTEST( test_complex(2*3*4*5*7) ); + CALL_SUBTEST( test_complex(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(256) ); + + CALL_SUBTEST( test_complex(3*8) ); + 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(5*32) ); + + CALL_SUBTEST( test_complex(2*3*4) ); + 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) ); + + CALL_SUBTEST( test_complex(2*3*4*5*7) ); + 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(32) ); - CALL_SUBTEST( test_scalar(45) ); CALL_SUBTEST( test_scalar(45) ); CALL_SUBTEST( test_scalar(45) ); - CALL_SUBTEST( test_scalar(50) ); CALL_SUBTEST( test_scalar(50) ); CALL_SUBTEST( test_scalar(50) ); - CALL_SUBTEST( test_scalar(256) ); 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) ); CALL_SUBTEST( test_scalar(2*3*4*5*7) ); + CALL_SUBTEST( test_scalar(32) ); + CALL_SUBTEST( test_scalar(32) ); + CALL_SUBTEST( test_scalar(32) ); + + CALL_SUBTEST( test_scalar(45) ); + CALL_SUBTEST( test_scalar(45) ); + CALL_SUBTEST( test_scalar(45) ); + + CALL_SUBTEST( test_scalar(50) ); + CALL_SUBTEST( test_scalar(50) ); + CALL_SUBTEST( test_scalar(50) ); + + CALL_SUBTEST( test_scalar(256) ); + 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) ); + CALL_SUBTEST( test_scalar(2*3*4*5*7) ); }