Update Eigen to commit:171bd08ca987987c3c50f0fa5dd8914bdd42dd3b
CHANGELOG
=========
171bd08ca - fix 2849
db85838ee - Add DUCC FFT support
6f1a14341 - Ensure info() implementation across all SolverBase derived types
f3e7d64f3 - Fix: Correct Lapacke bindings for BDCSVD and JacobiSVD to match the updated API
434a2fc4a - Fix obsolete comment in InverseImpl.h. We use PartialPivLU for the general case.
PiperOrigin-RevId: 759257247
Change-Id: I31deb8d802b88861782b34a089eecbdf71abfd47
diff --git a/Eigen/src/Core/SolverBase.h b/Eigen/src/Core/SolverBase.h
index df2ac83..5a6dfd4 100644
--- a/Eigen/src/Core/SolverBase.h
+++ b/Eigen/src/Core/SolverBase.h
@@ -78,6 +78,14 @@
template <typename Derived_>
friend struct internal::solve_assertion;
+ ComputationInfo info() const {
+ // CRTP static dispatch: Calls the 'info()' method on the derived class.
+ // Derived must implement 'ComputationInfo info() const'.
+ // If not implemented, name lookup falls back to this base method, causing
+ // infinite recursion (detectable by -Winfinite-recursion).
+ return derived().info();
+ }
+
enum {
RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
diff --git a/Eigen/src/Core/VectorwiseOp.h b/Eigen/src/Core/VectorwiseOp.h
index ac52dc5..9ccbf7d 100644
--- a/Eigen/src/Core/VectorwiseOp.h
+++ b/Eigen/src/Core/VectorwiseOp.h
@@ -603,10 +603,9 @@
/** Returns the expression where each subvector is the product of the vector \a other
* by the corresponding subvector of \c *this */
template <typename OtherDerived>
- EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
- CwiseBinaryOp<internal::scalar_product_op<Scalar>, const ExpressionTypeNestedCleaned,
- const typename ExtendedType<OtherDerived>::Type> EIGEN_DEVICE_FUNC
- operator*(const DenseBase<OtherDerived>& other) const {
+ EIGEN_DEVICE_FUNC CwiseBinaryOp<internal::scalar_product_op<Scalar, typename OtherDerived::Scalar>,
+ const ExpressionTypeNestedCleaned, const typename ExtendedType<OtherDerived>::Type>
+ operator*(const DenseBase<OtherDerived>& other) const {
EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType)
EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
@@ -616,8 +615,8 @@
/** Returns the expression where each subvector is the quotient of the corresponding
* subvector of \c *this by the vector \a other */
template <typename OtherDerived>
- EIGEN_DEVICE_FUNC CwiseBinaryOp<internal::scalar_quotient_op<Scalar>, const ExpressionTypeNestedCleaned,
- const typename ExtendedType<OtherDerived>::Type>
+ EIGEN_DEVICE_FUNC CwiseBinaryOp<internal::scalar_quotient_op<Scalar, typename OtherDerived::Scalar>,
+ const ExpressionTypeNestedCleaned, const typename ExtendedType<OtherDerived>::Type>
operator/(const DenseBase<OtherDerived>& other) const {
EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType)
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index a725a7b..786cd76 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -78,6 +78,17 @@
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime, PermutationIndex> PermutationPType;
typedef typename MatrixType::PlainObject PlainObject;
+ /** \brief Reports whether the LU factorization was successful.
+ *
+ * \note This function always returns \c Success. It is provided for compatibility
+ * with other factorization routines.
+ * \returns \c Success
+ */
+ ComputationInfo info() const {
+ eigen_assert(m_isInitialized && "FullPivLU is not initialized.");
+ return Success;
+ }
+
/**
* \brief Default Constructor.
*
diff --git a/Eigen/src/LU/InverseImpl.h b/Eigen/src/LU/InverseImpl.h
index 57fd677..fe8859e 100644
--- a/Eigen/src/LU/InverseImpl.h
+++ b/Eigen/src/LU/InverseImpl.h
@@ -268,7 +268,7 @@
* \note This matrix must be invertible, otherwise the result is undefined. If you need an
* invertibility check, do the following:
* \li for fixed sizes up to 4x4, use computeInverseAndDetWithCheck().
- * \li for the general case, use class FullPivLU.
+ * \li for the general case, use class PartialPivLU.
*
* Example: \include MatrixBase_inverse.cpp
* Output: \verbinclude MatrixBase_inverse.out
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index f09b90e..7ea14f5 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -90,6 +90,17 @@
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime, PermutationIndex> TranspositionType;
typedef typename MatrixType::PlainObject PlainObject;
+ /** \brief Reports whether the LU factorization was successful.
+ *
+ * \note This function always returns \c Success. It is provided for compatibility
+ * with other factorization routines.
+ * \returns \c Success
+ */
+ ComputationInfo info() const {
+ eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
+ return Success;
+ }
+
/**
* \brief Default Constructor.
*
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index cae9ae4..d173444 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -82,6 +82,17 @@
typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
typedef typename MatrixType::PlainObject PlainObject;
+ /** \brief Reports whether the QR factorization was successful.
+ *
+ * \note This function always returns \c Success. It is provided for compatibility
+ * with other factorization routines.
+ * \returns \c Success
+ */
+ ComputationInfo info() const {
+ eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
+ return Success;
+ }
+
/** \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index e297372..497085d 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -75,6 +75,17 @@
typedef HouseholderSequence<MatrixType, internal::remove_all_t<typename HCoeffsType::ConjugateReturnType>>
HouseholderSequenceType;
+ /** \brief Reports whether the QR factorization was successful.
+ *
+ * \note This function always returns \c Success. It is provided for compatibility
+ * with other factorization routines.
+ * \returns \c Success
+ */
+ ComputationInfo info() const {
+ eigen_assert(m_isInitialized && "HouseHolderQR is not initialized.");
+ return Success;
+ }
+
/**
* \brief Default Constructor.
*
diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h
index 85661a3..6fab905 100644
--- a/Eigen/src/SVD/BDCSVD.h
+++ b/Eigen/src/SVD/BDCSVD.h
@@ -182,7 +182,9 @@
* \deprecated Will be removed in the next major Eigen version. Options should
* be specified in the \a Options template parameter.
*/
- EIGEN_DEPRECATED BDCSVD(const MatrixType& matrix, unsigned int computationOptions) : m_algoswap(16), m_numIters(0) {
+ template <typename Derived>
+ EIGEN_DEPRECATED BDCSVD(const MatrixBase<Derived>& matrix, unsigned int computationOptions)
+ : m_algoswap(16), m_numIters(0) {
internal::check_svd_options_assertions<MatrixType, Options>(computationOptions, matrix.rows(), matrix.cols());
compute_impl(matrix, computationOptions);
}
diff --git a/Eigen/src/SVD/BDCSVD_LAPACKE.h b/Eigen/src/SVD/BDCSVD_LAPACKE.h
index 89d5cbd..5d2b8c7 100644
--- a/Eigen/src/SVD/BDCSVD_LAPACKE.h
+++ b/Eigen/src/SVD/BDCSVD_LAPACKE.h
@@ -58,7 +58,8 @@
// construct this by moving from a parent object
BDCSVD_LAPACKE(SVD&& svd) : SVD(std::move(svd)) {}
- void compute_impl_lapacke(const MatrixType& matrix, unsigned int computationOptions) {
+ template <typename Derived>
+ void compute_impl_lapacke(const MatrixBase<Derived>& matrix, unsigned int computationOptions) {
SVD::allocate(matrix.rows(), matrix.cols(), computationOptions);
SVD::m_nonzeroSingularValues = SVD::m_diagSize;
@@ -120,8 +121,8 @@
}
};
-template <typename MatrixType_, int Options>
-BDCSVD<MatrixType_, Options>& BDCSVD_wrapper(BDCSVD<MatrixType_, Options>& svd, const MatrixType_& matrix,
+template <typename MatrixType_, int Options, typename Derived>
+BDCSVD<MatrixType_, Options>& BDCSVD_wrapper(BDCSVD<MatrixType_, Options>& svd, const MatrixBase<Derived>& matrix,
int computationOptions) {
// we need to move to the wrapper type and back
BDCSVD_LAPACKE<MatrixType_, Options> tmpSvd(std::move(svd));
@@ -134,12 +135,13 @@
} // end namespace internal
-#define EIGEN_LAPACKE_SDD(EIGTYPE, EIGCOLROW, OPTIONS) \
- template <> \
- inline BDCSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>& \
- BDCSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>::compute_impl( \
- const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>& matrix, unsigned int computationOptions) { \
- return internal::lapacke_helpers::BDCSVD_wrapper(*this, matrix, computationOptions); \
+#define EIGEN_LAPACKE_SDD(EIGTYPE, EIGCOLROW, OPTIONS) \
+ template <> \
+ template <typename Derived> \
+ inline BDCSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>& \
+ BDCSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>::compute_impl( \
+ const MatrixBase<Derived>& matrix, unsigned int computationOptions) { \
+ return internal::lapacke_helpers::BDCSVD_wrapper(*this, matrix, computationOptions); \
}
#define EIGEN_LAPACK_SDD_OPTIONS(OPTIONS) \
diff --git a/Eigen/src/SVD/JacobiSVD_LAPACKE.h b/Eigen/src/SVD/JacobiSVD_LAPACKE.h
index df6a096..db26366 100644
--- a/Eigen/src/SVD/JacobiSVD_LAPACKE.h
+++ b/Eigen/src/SVD/JacobiSVD_LAPACKE.h
@@ -40,65 +40,65 @@
/** \internal Specialization for the data types supported by LAPACKe */
-#define EIGEN_LAPACKE_SVD(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW, OPTIONS) \
- template <> \
- inline JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>& \
- JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>::compute_impl( \
- const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>& matrix, unsigned int computationOptions) { \
- typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> MatrixType; \
- /*typedef MatrixType::Scalar Scalar;*/ \
- /*typedef MatrixType::RealScalar RealScalar;*/ \
- allocate(matrix.rows(), matrix.cols(), computationOptions); \
- \
- /*const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();*/ \
- m_nonzeroSingularValues = diagSize(); \
- \
- lapack_int lda = internal::convert_index<lapack_int>(matrix.outerStride()), ldu, ldvt; \
- lapack_int matrix_order = LAPACKE_COLROW; \
- char jobu, jobvt; \
- LAPACKE_TYPE *u, *vt, dummy; \
- jobu = (m_computeFullU) ? 'A' : (m_computeThinU) ? 'S' : 'N'; \
- jobvt = (m_computeFullV) ? 'A' : (m_computeThinV) ? 'S' : 'N'; \
- if (computeU()) { \
- ldu = internal::convert_index<lapack_int>(m_matrixU.outerStride()); \
- u = (LAPACKE_TYPE*)m_matrixU.data(); \
- } else { \
- ldu = 1; \
- u = &dummy; \
- } \
- MatrixType localV; \
- lapack_int vt_rows = (m_computeFullV) ? internal::convert_index<lapack_int>(cols()) \
- : (m_computeThinV) ? internal::convert_index<lapack_int>(diagSize()) \
- : 1; \
- if (computeV()) { \
- localV.resize(vt_rows, cols()); \
- ldvt = internal::convert_index<lapack_int>(localV.outerStride()); \
- vt = (LAPACKE_TYPE*)localV.data(); \
- } else { \
- ldvt = 1; \
- vt = &dummy; \
- } \
- Matrix<LAPACKE_RTYPE, Dynamic, Dynamic> superb; \
- superb.resize(diagSize(), 1); \
- MatrixType m_temp; \
- m_temp = matrix; \
- lapack_int info = LAPACKE_##LAPACKE_PREFIX##gesvd( \
- matrix_order, jobu, jobvt, internal::convert_index<lapack_int>(rows()), \
- internal::convert_index<lapack_int>(cols()), (LAPACKE_TYPE*)m_temp.data(), lda, \
- (LAPACKE_RTYPE*)m_singularValues.data(), u, ldu, vt, ldvt, superb.data()); \
- /* Check the result of the LAPACK call */ \
- if (info < 0 || !m_singularValues.allFinite()) { \
- m_info = InvalidInput; \
- } else if (info > 0) { \
- m_info = NoConvergence; \
- } else { \
- m_info = Success; \
- if (computeV()) m_matrixV = localV.adjoint(); \
- } \
- /* for(int i=0;i<diagSize();i++) if (m_singularValues.coeffRef(i) < precision) { m_nonzeroSingularValues--; \
- * m_singularValues.coeffRef(i)=RealScalar(0);}*/ \
- m_isInitialized = true; \
- return *this; \
+#define EIGEN_LAPACKE_SVD(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW, OPTIONS) \
+ template <> \
+ template <typename Derived> \
+ inline JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>& \
+ JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>::compute_impl( \
+ const MatrixBase<Derived>& matrix, unsigned int computationOptions) { \
+ /*typedef MatrixType::Scalar Scalar;*/ \
+ /*typedef MatrixType::RealScalar RealScalar;*/ \
+ allocate(matrix.rows(), matrix.cols(), computationOptions); \
+ \
+ /*const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();*/ \
+ m_nonzeroSingularValues = diagSize(); \
+ \
+ lapack_int lda = internal::convert_index<lapack_int>(matrix.outerStride()), ldu, ldvt; \
+ lapack_int matrix_order = LAPACKE_COLROW; \
+ char jobu, jobvt; \
+ LAPACKE_TYPE *u, *vt, dummy; \
+ jobu = (m_computeFullU) ? 'A' : (m_computeThinU) ? 'S' : 'N'; \
+ jobvt = (m_computeFullV) ? 'A' : (m_computeThinV) ? 'S' : 'N'; \
+ if (computeU()) { \
+ ldu = internal::convert_index<lapack_int>(m_matrixU.outerStride()); \
+ u = (LAPACKE_TYPE*)m_matrixU.data(); \
+ } else { \
+ ldu = 1; \
+ u = &dummy; \
+ } \
+ MatrixType localV; \
+ lapack_int vt_rows = (m_computeFullV) ? internal::convert_index<lapack_int>(cols()) \
+ : (m_computeThinV) ? internal::convert_index<lapack_int>(diagSize()) \
+ : 1; \
+ if (computeV()) { \
+ localV.resize(vt_rows, cols()); \
+ ldvt = internal::convert_index<lapack_int>(localV.outerStride()); \
+ vt = (LAPACKE_TYPE*)localV.data(); \
+ } else { \
+ ldvt = 1; \
+ vt = &dummy; \
+ } \
+ Matrix<LAPACKE_RTYPE, Dynamic, Dynamic> superb; \
+ superb.resize(diagSize(), 1); \
+ MatrixType m_temp; \
+ m_temp = matrix; \
+ lapack_int info = LAPACKE_##LAPACKE_PREFIX##gesvd( \
+ matrix_order, jobu, jobvt, internal::convert_index<lapack_int>(rows()), \
+ internal::convert_index<lapack_int>(cols()), (LAPACKE_TYPE*)m_temp.data(), lda, \
+ (LAPACKE_RTYPE*)m_singularValues.data(), u, ldu, vt, ldvt, superb.data()); \
+ /* Check the result of the LAPACK call */ \
+ if (info < 0 || !m_singularValues.allFinite()) { \
+ m_info = InvalidInput; \
+ } else if (info > 0) { \
+ m_info = NoConvergence; \
+ } else { \
+ m_info = Success; \
+ if (computeV()) m_matrixV = localV.adjoint(); \
+ } \
+ /* for(int i=0;i<diagSize();i++) if (m_singularValues.coeffRef(i) < precision) { m_nonzeroSingularValues--; \
+ * m_singularValues.coeffRef(i)=RealScalar(0);}*/ \
+ m_isInitialized = true; \
+ return *this; \
}
#define EIGEN_LAPACK_SVD_OPTIONS(OPTIONS) \
diff --git a/test/vectorwiseop.cpp b/test/vectorwiseop.cpp
index 0d1896b..6d0e5cb 100644
--- a/test/vectorwiseop.cpp
+++ b/test/vectorwiseop.cpp
@@ -214,6 +214,17 @@
VERIFY_IS_EQUAL(m1.real().middleCols(0, fix<0>).colwise().maxCoeff().eval().cols(), 0);
}
+void vectorwiseop_mixedscalar() {
+ Matrix4cd a = Matrix4cd::Random();
+ Vector4cd b = Vector4cd::Random();
+ b.imag().setZero();
+ Vector4d b_real = b.real();
+
+ Matrix4cd c = a.array().rowwise() * b.array().transpose();
+ Matrix4cd d = a.array().rowwise() * b_real.array().transpose();
+ VERIFY_IS_CWISE_EQUAL(c, d);
+}
+
EIGEN_DECLARE_TEST(vectorwiseop) {
CALL_SUBTEST_1(vectorwiseop_array(Array22cd()));
CALL_SUBTEST_2(vectorwiseop_array(Array<double, 3, 2>()));
@@ -226,4 +237,5 @@
MatrixXd(internal::random<int>(1, EIGEN_TEST_MAX_SIZE), internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
CALL_SUBTEST_7(vectorwiseop_matrix(VectorXd(internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
CALL_SUBTEST_7(vectorwiseop_matrix(RowVectorXd(internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
+ CALL_SUBTEST_8(vectorwiseop_mixedscalar());
}
diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT
index f32ad36..bca5a3e 100644
--- a/unsupported/Eigen/FFT
+++ b/unsupported/Eigen/FFT
@@ -37,12 +37,13 @@
* - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
* - MKL (https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-download.html) : fastest, free -- may be
* incompatible with Eigen in GPL form.
- * - pocketfft (https://gitlab.mpcdf.mpg.de/mtr/pocketfft) : faster than kissfft, BSD 3-clause.
+ * - PocketFFT/DUCC (https://gitlab.mpcdf.mpg.de/mtr/pocketfft, https://gitlab.mpcdf.mpg.de/mtr/ducc) : 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
+ * According to the author, DUCC contains the "evolution" of pocketfft, though the interface is very similar.
* used for these cases
*
* \section FFTDesign Design
@@ -85,7 +86,7 @@
#ifdef EIGEN_FFTW_DEFAULT
// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
#include <fftw3.h>
-#include "src/FFT/ei_fftw_impl.h"
+#include "src/FFT/fftw_impl.h"
namespace Eigen {
// template <typename T> typedef struct internal::fftw_impl default_fft_impl; this does not work
template <typename T>
@@ -93,7 +94,7 @@
} // namespace Eigen
#elif defined EIGEN_MKL_DEFAULT
// intel Math Kernel Library: fastest, free -- may be incompatible with Eigen in GPL form
-#include "src/FFT/ei_imklfft_impl.h"
+#include "src/FFT/imklfft_impl.h"
namespace Eigen {
template <typename T>
struct default_fft_impl : public internal::imklfft::imklfft_impl<T> {};
@@ -101,14 +102,24 @@
#elif defined EIGEN_POCKETFFT_DEFAULT
// internal::pocketfft_impl: a heavily modified implementation of FFTPack, with many advantages.
#include <pocketfft_hdronly.h>
-#include "src/FFT/ei_pocketfft_impl.h"
+#include "src/FFT/pocketfft_impl.h"
namespace Eigen {
template <typename T>
struct default_fft_impl : public internal::pocketfft_impl<T> {};
} // namespace Eigen
+#elif defined EIGEN_DUCCFFT_DEFAULT
+#include <ducc0/fft/fft.h>
+#include <ducc0/infra/string_utils.h>
+#include <ducc0/fft/fft.h>
+#include <ducc0/fft/fftnd_impl.h>
+#include "src/FFT/duccfft_impl.h"
+namespace Eigen {
+template <typename T>
+struct default_fft_impl : public internal::duccfft_impl<T> {};
+} // namespace Eigen
#else
// internal::kissfft_impl: small, free, reasonably efficient default, derived from kissfft
-#include "src/FFT/ei_kissfft_impl.h"
+#include "src/FFT/kissfft_impl.h"
namespace Eigen {
template <typename T>
struct default_fft_impl : public internal::kissfft_impl<T> {};
@@ -204,7 +215,8 @@
inline void fwd(Complex* dst, const Complex* src, Index nfft) { m_impl.fwd(dst, src, static_cast<int>(nfft)); }
-#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
+#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT || \
+ defined EIGEN_MKL_DEFAULT
inline void fwd2(Complex* dst, const Complex* src, int n0, int n1) { m_impl.fwd2(dst, src, n0, n1); }
#endif
@@ -366,7 +378,8 @@
inv(&dst[0], &src[0], nfft);
}
-#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
+#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT || \
+ defined EIGEN_MKL_DEFAULT
inline void inv2(Complex* dst, const Complex* src, int n0, int n1) {
m_impl.inv2(dst, src, n0, n1);
if (HasFlag(Unscaled) == false) scale(dst, 1. / (n0 * n1), n0 * n1);
@@ -385,7 +398,6 @@
Matrix<T_Data, Dynamic, 1>::Map(x, nx) *= s;
else
Matrix<T_Data, Dynamic, 1>::MapAligned(x, nx) *= s;
- // Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s;
#endif
}
diff --git a/unsupported/Eigen/src/FFT/duccfft_impl.h b/unsupported/Eigen/src/FFT/duccfft_impl.h
new file mode 100644
index 0000000..781716d
--- /dev/null
+++ b/unsupported/Eigen/src/FFT/duccfft_impl.h
@@ -0,0 +1,71 @@
+// 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/.
+
+namespace Eigen {
+
+namespace internal {
+
+template <typename _Scalar>
+struct duccfft_impl {
+ using Scalar = _Scalar;
+ using Complex = std::complex<Scalar>;
+ using shape_t = ducc0::fmav_info::shape_t;
+ using stride_t = ducc0::fmav_info::stride_t;
+
+ inline void clear() {}
+
+ inline void fwd(Complex* dst, const Scalar* src, int nfft) {
+ const shape_t axes{0};
+ ducc0::cfmav<Scalar> m_in(src, shape_t{static_cast<size_t>(nfft)});
+ ducc0::vfmav<Complex> m_out(dst, shape_t{static_cast<size_t>(nfft) / 2 + 1});
+ ducc0::r2c(m_in, m_out, axes, /*forward=*/true, /*scale=*/static_cast<Scalar>(1));
+ }
+
+ inline void fwd(Complex* dst, const Complex* src, int nfft) {
+ const shape_t axes{0};
+ ducc0::cfmav<Complex> m_in(src, shape_t{static_cast<size_t>(nfft)});
+ ducc0::vfmav<Complex> m_out(dst, shape_t{static_cast<size_t>(nfft)});
+ ducc0::c2c(m_in, m_out, axes, /*forward=*/true, /*scale=*/static_cast<Scalar>(1));
+ }
+
+ inline void inv(Scalar* dst, const Complex* src, int nfft) {
+ const shape_t axes{0};
+ ducc0::cfmav<Complex> m_in(src, shape_t{static_cast<size_t>(nfft) / 2 + 1});
+ ducc0::vfmav<Scalar> m_out(dst, shape_t{static_cast<size_t>(nfft)});
+ ducc0::c2r(m_in, m_out, axes, /*forward=*/false, /*scale=*/static_cast<Scalar>(1));
+ }
+
+ inline void inv(Complex* dst, const Complex* src, int nfft) {
+ const shape_t axes{0};
+ ducc0::cfmav<Complex> m_in(src, shape_t{static_cast<size_t>(nfft)});
+ ducc0::vfmav<Complex> m_out(dst, shape_t{static_cast<size_t>(nfft)});
+ ducc0::c2c(m_in, m_out, axes, /*forward=*/false, /*scale=*/static_cast<Scalar>(1));
+ }
+
+ inline void fwd2(Complex* dst, const Complex* src, int nfft0, int nfft1) {
+ const shape_t axes{0, 1};
+ const shape_t in_shape{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
+ const shape_t out_shape{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
+ const stride_t stride{static_cast<ptrdiff_t>(nfft1), static_cast<ptrdiff_t>(1)};
+ ducc0::cfmav<Complex> m_in(src, in_shape, stride);
+ ducc0::vfmav<Complex> m_out(dst, out_shape, stride);
+ ducc0::c2c(m_in, m_out, axes, /*forward=*/true, /*scale=*/static_cast<Scalar>(1));
+ }
+
+ inline void inv2(Complex* dst, const Complex* src, int nfft0, int nfft1) {
+ const shape_t axes{0, 1};
+ const shape_t in_shape{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
+ const shape_t out_shape{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
+ const stride_t stride{static_cast<ptrdiff_t>(nfft1), static_cast<ptrdiff_t>(1)};
+ ducc0::cfmav<Complex> m_in(src, in_shape, stride);
+ ducc0::vfmav<Complex> m_out(dst, out_shape, stride);
+ ducc0::c2c(m_in, m_out, axes, /*forward=*/false, /*scale=*/static_cast<Scalar>(1));
+ }
+};
+
+} // namespace internal
+} // namespace Eigen
diff --git a/unsupported/Eigen/src/FFT/ei_fftw_impl.h b/unsupported/Eigen/src/FFT/fftw_impl.h
similarity index 100%
rename from unsupported/Eigen/src/FFT/ei_fftw_impl.h
rename to unsupported/Eigen/src/FFT/fftw_impl.h
diff --git a/unsupported/Eigen/src/FFT/ei_imklfft_impl.h b/unsupported/Eigen/src/FFT/imklfft_impl.h
similarity index 100%
rename from unsupported/Eigen/src/FFT/ei_imklfft_impl.h
rename to unsupported/Eigen/src/FFT/imklfft_impl.h
diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/kissfft_impl.h
similarity index 100%
rename from unsupported/Eigen/src/FFT/ei_kissfft_impl.h
rename to unsupported/Eigen/src/FFT/kissfft_impl.h
diff --git a/unsupported/Eigen/src/FFT/ei_pocketfft_impl.h b/unsupported/Eigen/src/FFT/pocketfft_impl.h
similarity index 71%
rename from unsupported/Eigen/src/FFT/ei_pocketfft_impl.h
rename to unsupported/Eigen/src/FFT/pocketfft_impl.h
index 918ceb0..fce105b 100644
--- a/unsupported/Eigen/src/FFT/ei_pocketfft_impl.h
+++ b/unsupported/Eigen/src/FFT/pocketfft_impl.h
@@ -5,17 +5,16 @@
// 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 <typename _Scalar>
struct pocketfft_impl {
- typedef _Scalar Scalar;
- typedef std::complex<Scalar> Complex;
+ using Scalar = _Scalar;
+ using Complex = std::complex<Scalar>;
+ using shape_t = pocketfft::shape_t;
+ using stride_t = pocketfft::stride_t;
inline void clear() {}
@@ -24,14 +23,14 @@
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<Scalar>(1));
+ pocketfft::r2c(shape_, stride_in, stride_out, axes_, pocketfft::FORWARD, src, dst, static_cast<Scalar>(1));
}
inline void fwd(Complex* dst, const Complex* src, int nfft) {
const shape_t shape_{static_cast<size_t>(nfft)};
const shape_t axes_{0};
const stride_t stride_{sizeof(Complex)};
- c2c(shape_, stride_, stride_, axes_, FORWARD, src, dst, static_cast<Scalar>(1));
+ pocketfft::c2c(shape_, stride_, stride_, axes_, pocketfft::FORWARD, src, dst, static_cast<Scalar>(1));
}
inline void inv(Scalar* dst, const Complex* src, int nfft) {
@@ -39,28 +38,28 @@
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<Scalar>(1));
+ pocketfft::c2r(shape_, stride_in, stride_out, axes_, pocketfft::BACKWARD, src, dst, static_cast<Scalar>(1));
}
inline void inv(Complex* dst, const Complex* src, int nfft) {
const shape_t shape_{static_cast<size_t>(nfft)};
const shape_t axes_{0};
const stride_t stride_{sizeof(Complex)};
- c2c(shape_, stride_, stride_, axes_, BACKWARD, src, dst, static_cast<Scalar>(1));
+ pocketfft::c2c(shape_, stride_, stride_, axes_, pocketfft::BACKWARD, src, dst, static_cast<Scalar>(1));
}
inline void fwd2(Complex* dst, const Complex* src, int nfft0, int nfft1) {
const shape_t shape_{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
const shape_t axes_{0, 1};
const stride_t stride_{static_cast<ptrdiff_t>(sizeof(Complex) * nfft1), static_cast<ptrdiff_t>(sizeof(Complex))};
- c2c(shape_, stride_, stride_, axes_, FORWARD, src, dst, static_cast<Scalar>(1));
+ pocketfft::c2c(shape_, stride_, stride_, axes_, pocketfft::FORWARD, src, dst, static_cast<Scalar>(1));
}
inline void inv2(Complex* dst, const Complex* src, int nfft0, int nfft1) {
const shape_t shape_{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
const shape_t axes_{0, 1};
const stride_t stride_{static_cast<ptrdiff_t>(sizeof(Complex) * nfft1), static_cast<ptrdiff_t>(sizeof(Complex))};
- c2c(shape_, stride_, stride_, axes_, BACKWARD, src, dst, static_cast<Scalar>(1));
+ pocketfft::c2c(shape_, stride_, stride_, axes_, pocketfft::BACKWARD, src, dst, static_cast<Scalar>(1));
}
};
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index c270458..5efa7e8 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -88,6 +88,25 @@
ei_add_property(EIGEN_MISSING_BACKENDS "pocketfft, ")
endif()
+if( NOT DUCC_ROOT AND ENV{DUCC_ROOT} )
+ set( DUCC_ROOT $ENV{DUCC_ROOT} )
+endif()
+find_path(DUCCFFT
+ NAMES "src/ducc0/fft/fft.h"
+ PATHS ${DUCC_ROOT})
+message(INFO " ${DUCC_ROOT} ${DUCCFFT}")
+if(DUCCFFT)
+ ei_add_property(EIGEN_TESTED_BACKENDS "duccfft, ")
+ include_directories( "${DUCCFFT}/src" )
+ add_library(ducc_lib "${DUCCFFT}/src/ducc0/infra/string_utils.cc" "${DUCCFFT}/src/ducc0/infra/threading.cc")
+ target_compile_definitions(ducc_lib PUBLIC "DUCC0_NO_THREADING=1")
+ ei_add_test(duccfft "-DEIGEN_DUCCFFT_DEFAULT -DDUCC0_NO_THREADING=1" "ducc_lib" )
+ set_target_properties(ducc_lib duccfft PROPERTIES CXX_STANDARD 17)
+else()
+ ei_add_property(EIGEN_MISSING_BACKENDS "duccfft, ")
+endif()
+
+
option(EIGEN_TEST_OPENGL "Enable OpenGL support in unit tests" OFF)
if(EIGEN_TEST_OPENGL)
find_package(OpenGL)
diff --git a/unsupported/test/duccfft.cpp b/unsupported/test/duccfft.cpp
new file mode 100644
index 0000000..4247663
--- /dev/null
+++ b/unsupported/test/duccfft.cpp
@@ -0,0 +1,4 @@
+#define EIGEN_DUCCFFT_DEFAULT 1
+#include <ducc0/fft/fft.h> // Needs to be included before main.h
+#include <ducc0/fft/fftnd_impl.h> // Same requirement
+#include "fft_test_shared.h"
diff --git a/unsupported/test/fft_test_shared.h b/unsupported/test/fft_test_shared.h
index 3adcd90..d7380be 100644
--- a/unsupported/test/fft_test_shared.h
+++ b/unsupported/test/fft_test_shared.h
@@ -272,7 +272,7 @@
CALL_SUBTEST(test_scalar<float>(2 * 3 * 4 * 5 * 7));
CALL_SUBTEST(test_scalar<double>(2 * 3 * 4 * 5 * 7));
-#if defined EIGEN_HAS_FFTWL || defined EIGEN_POCKETFFT_DEFAULT
+#if defined EIGEN_HAS_FFTWL || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT
CALL_SUBTEST(test_complex<long double>(32));
CALL_SUBTEST(test_complex<long double>(256));
CALL_SUBTEST(test_complex<long double>(3 * 8));
@@ -294,13 +294,15 @@
// fail to build since Eigen limit the stack allocation size,too big here.
// CALL_SUBTEST( ( test_complex2d<long double, 256, 256> () ) );
#endif
-#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
+#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT || \
+ defined EIGEN_MKL_DEFAULT
CALL_SUBTEST((test_complex2d<float, 24, 24>()));
CALL_SUBTEST((test_complex2d<float, 60, 60>()));
CALL_SUBTEST((test_complex2d<float, 24, 60>()));
CALL_SUBTEST((test_complex2d<float, 60, 24>()));
#endif
-#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
+#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT || \
+ defined EIGEN_MKL_DEFAULT
CALL_SUBTEST((test_complex2d<double, 24, 24>()));
CALL_SUBTEST((test_complex2d<double, 60, 60>()));
CALL_SUBTEST((test_complex2d<double, 24, 60>()));
