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>()));