Core: Vectorize complex component reductions libeigen/eigen!3108 Closes #3097 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h
index d4f351f..ff9a37f 100644
--- a/Eigen/src/Core/CoreEvaluators.h
+++ b/Eigen/src/Core/CoreEvaluators.h

@@ -1169,6 +1169,100 @@
   Data m_d;
 };
 
+// Read-only component evaluators can gather from contiguous std::complex storage. In particular,
+// PacketAccessBit on a custom expression does not imply array-oriented component access.
+template <typename XprType, int Component,
+          bool Vectorizable = bool(traits<typename XprType::NestedExpression>::Flags & DirectAccessBit) &&
+                              inner_stride_at_compile_time<typename XprType::NestedExpression>::value == 1 &&
+                              (std::is_same<typename XprType::Scalar, float>::value ||
+                               std::is_same<typename XprType::Scalar, double>::value) &&
+                              packet_traits<typename XprType::Scalar>::Vectorizable &&
+                              std::is_same<typename XprType::Scalar, typename unpacket_traits<typename packet_traits<
+                                                                         typename XprType::Scalar>::type>::type>::value>
+struct complex_component_evaluator : unary_evaluator<XprType> {
+  using Base = unary_evaluator<XprType>;
+  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE explicit complex_component_evaluator(const XprType& xpr)
+      : Base(xpr) {}
+};
+
+template <typename XprType, int Component>
+struct complex_component_evaluator<XprType, Component, true> : unary_evaluator<XprType> {
+  using Base = unary_evaluator<XprType>;
+  using Scalar = typename XprType::Scalar;
+  using ArgType = typename XprType::NestedExpression;
+  // Blocks must delegate packet reads here rather than reinterpret this stride-2 view as a contiguous Map.
+  static constexpr unsigned int Flags = (Base::Flags & ~DirectAccessBit) | PacketAccessBit;
+  static constexpr int Alignment = 0;
+
+  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE explicit complex_component_evaluator(const XprType& xpr)
+      : Base(xpr),
+        m_data(reinterpret_cast<const Scalar*>(xpr.nestedExpression().data())),
+        m_outerStride(xpr.nestedExpression().outerStride()) {}
+
+  template <int LoadMode, typename PacketType>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PacketType packet(Index index) const {
+    return pgather<Scalar, PacketType>(m_data + 2 * index + Component, 2);
+  }
+
+  template <int LoadMode, typename PacketType>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const {
+    const Index index = XprType::IsRowMajor ? row * m_outerStride.value() + col : col * m_outerStride.value() + row;
+    return packet<LoadMode, PacketType>(index);
+  }
+
+  template <int LoadMode, typename PacketType>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PacketType packetSegment(Index index, Index begin, Index count) const {
+    Scalar values[unpacket_traits<PacketType>::size] = {};
+    for (Index i = begin; i < begin + count; ++i) values[i] = Base::coeff(index + i);
+    return ploadu<PacketType>(values);
+  }
+
+  template <int LoadMode, typename PacketType>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PacketType packetSegment(Index row, Index col, Index begin, Index count) const {
+    Scalar values[unpacket_traits<PacketType>::size] = {};
+    for (Index i = begin; i < begin + count; ++i)
+      values[i] = Base::coeff(row + (XprType::IsRowMajor ? 0 : i), col + (XprType::IsRowMajor ? i : 0));
+    return ploadu<PacketType>(values);
+  }
+
+ private:
+  const Scalar* m_data;
+  const variable_if_dynamic<Index, outer_stride_at_compile_time<ArgType>::value> m_outerStride;
+};
+
+template <typename Real, typename ArgType>
+struct evaluator<CwiseUnaryOp<scalar_real_op<std::complex<Real>>, ArgType>>
+    : complex_component_evaluator<CwiseUnaryOp<scalar_real_op<std::complex<Real>>, ArgType>, 0> {
+  using XprType = CwiseUnaryOp<scalar_real_op<std::complex<Real>>, ArgType>;
+  using Base = complex_component_evaluator<XprType, 0>;
+  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr) : Base(xpr) {}
+};
+
+// Writable views keep scalar stores: assignment functors require contiguous destination packets.
+template <typename Real, typename ArgType>
+struct evaluator<const CwiseUnaryView<scalar_real_ref_op<std::complex<Real>>, ArgType, Stride<0, 0>>>
+    : complex_component_evaluator<CwiseUnaryView<scalar_real_ref_op<std::complex<Real>>, ArgType, Stride<0, 0>>, 0> {
+  using XprType = CwiseUnaryView<scalar_real_ref_op<std::complex<Real>>, ArgType, Stride<0, 0>>;
+  using Base = complex_component_evaluator<XprType, 0>;
+  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr) : Base(xpr) {}
+};
+
+template <typename Real, typename ArgType>
+struct evaluator<CwiseUnaryOp<scalar_imag_op<std::complex<Real>>, ArgType>>
+    : complex_component_evaluator<CwiseUnaryOp<scalar_imag_op<std::complex<Real>>, ArgType>, 1> {
+  using XprType = CwiseUnaryOp<scalar_imag_op<std::complex<Real>>, ArgType>;
+  using Base = complex_component_evaluator<XprType, 1>;
+  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr) : Base(xpr) {}
+};
+
+template <typename Real, typename ArgType>
+struct evaluator<const CwiseUnaryView<scalar_imag_ref_op<std::complex<Real>>, ArgType, Stride<0, 0>>>
+    : complex_component_evaluator<CwiseUnaryView<scalar_imag_ref_op<std::complex<Real>>, ArgType, Stride<0, 0>>, 1> {
+  using XprType = CwiseUnaryView<scalar_imag_ref_op<std::complex<Real>>, ArgType, Stride<0, 0>>;
+  using Base = complex_component_evaluator<XprType, 1>;
+  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr) : Base(xpr) {}
+};
+
 // -------------------- Map --------------------
 
 // FIXME: consider using Derived::PlainObject for PlainObjectType.
@@ -1325,7 +1419,7 @@
 // -------------------- Block --------------------
 
 template <typename ArgType, int BlockRows, int BlockCols, bool InnerPanel,
-          bool HasDirectAccess = has_direct_access<ArgType>::value>
+          bool HasDirectAccess = has_direct_access<ArgType>::value && bool(evaluator<ArgType>::Flags & DirectAccessBit)>
 struct block_evaluator;
 
 template <typename ArgType, int BlockRows, int BlockCols, bool InnerPanel>
@@ -1354,7 +1448,11 @@
                                                             : int(outer_stride_at_compile_time<ArgType>::value),
     OuterStrideAtCompileTime = HasSameStorageOrderAsArgType ? int(outer_stride_at_compile_time<ArgType>::value)
                                                             : int(inner_stride_at_compile_time<ArgType>::value),
-    MaskPacketAccessBit = (InnerStrideAtCompileTime == 1 || HasSameStorageOrderAsArgType) ? PacketAccessBit : 0,
+    // Direct-access blocks use mapbase_evaluator's contiguous loads, not the parent's packet method.
+    MaskPacketAccessBit = (InnerStrideAtCompileTime == 1 ||
+                           (!(evaluator<ArgType>::Flags & DirectAccessBit) && HasSameStorageOrderAsArgType))
+                              ? PacketAccessBit
+                              : 0,
 
     FlagsLinearAccessBit = (RowsAtCompileTime == 1 || ColsAtCompileTime == 1 ||
                             (InnerPanel && (evaluator<ArgType>::Flags & LinearAccessBit)))

diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h
index c19a00a..133cbdc 100644
--- a/Eigen/src/Core/Redux.h
+++ b/Eigen/src/Core/Redux.h

@@ -492,8 +492,8 @@
 
 // evaluator adaptor
 template <typename XprType_>
-class redux_evaluator : public internal::evaluator<XprType_> {
-  using Base = internal::evaluator<XprType_>;
+class redux_evaluator : public internal::evaluator<const XprType_> {
+  using Base = internal::evaluator<const XprType_>;
 
  public:
   using XprType = XprType_;

diff --git a/Eigen/src/Core/arch/AVX/PacketMath.h b/Eigen/src/Core/arch/AVX/PacketMath.h
index e7a4dc4..129c8e0 100644
--- a/Eigen/src/Core/arch/AVX/PacketMath.h
+++ b/Eigen/src/Core/arch/AVX/PacketMath.h

@@ -1638,11 +1638,32 @@
 // 4);
 template <>
 EIGEN_DEVICE_FUNC inline Packet8f pgather<float, Packet8f>(const float* from, Index stride) {
+  if (stride == 2) {
+#ifdef EIGEN_VECTORIZE_AVX2
+    // Both loads stay within from[0..14], including when from points to an imaginary component.
+    const Packet8f low = _mm256_loadu_ps(from);
+    const Packet8f high = _mm256_loadu_ps(from + 7);
+    const Packet8f interleaved = _mm256_shuffle_ps(low, high, _MM_SHUFFLE(3, 1, 2, 0));
+    return _mm256_castpd_ps(_mm256_permute4x64_pd(_mm256_castps_pd(interleaved), _MM_SHUFFLE(3, 1, 2, 0)));
+#else
+    const Packet4f low = pgather<float, Packet4f>(from, 2);
+    const Packet4f high = pgather<float, Packet4f>(from + 8, 2);
+    return _mm256_insertf128_ps(_mm256_castps128_ps256(low), high, 1);
+#endif
+  }
   return _mm256_set_ps(from[7 * stride], from[6 * stride], from[5 * stride], from[4 * stride], from[3 * stride],
                        from[2 * stride], from[1 * stride], from[0 * stride]);
 }
 template <>
 EIGEN_DEVICE_FUNC inline Packet4d pgather<double, Packet4d>(const double* from, Index stride) {
+#ifdef EIGEN_VECTORIZE_AVX2
+  if (stride == 2) {
+    const Packet4d low = _mm256_loadu_pd(from);
+    const Packet4d high = _mm256_loadu_pd(from + 3);
+    const Packet4d interleaved = _mm256_shuffle_pd(low, high, 0xa);
+    return _mm256_permute4x64_pd(interleaved, _MM_SHUFFLE(3, 1, 2, 0));
+  }
+#endif
   return _mm256_set_pd(from[3 * stride], from[2 * stride], from[1 * stride], from[0 * stride]);
 }
 template <>

diff --git a/Eigen/src/Core/arch/AVX512/PacketMath.h b/Eigen/src/Core/arch/AVX512/PacketMath.h
index 3f1ba2d..2596313 100644
--- a/Eigen/src/Core/arch/AVX512/PacketMath.h
+++ b/Eigen/src/Core/arch/AVX512/PacketMath.h

@@ -1247,6 +1247,13 @@
 
 template <>
 EIGEN_DEVICE_FUNC inline Packet16f pgather<float, Packet16f>(const float* from, Index stride) {
+  if (stride == 2) {
+    // The overlapping loads end at from[30], not the unused coefficient from[31].
+    const Packet16f low = _mm512_loadu_ps(from);
+    const Packet16f high = _mm512_loadu_ps(from + 15);
+    const Packet16i indices = _mm512_setr_epi32(0, 2, 4, 6, 8, 10, 12, 14, 17, 19, 21, 23, 25, 27, 29, 31);
+    return _mm512_permutex2var_ps(low, indices, high);
+  }
   Packet16i stride_vector = _mm512_set1_epi32(convert_index<int>(stride));
   Packet16i stride_multiplier = _mm512_set_epi32(15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0);
   Packet16i indices = _mm512_mullo_epi32(stride_vector, stride_multiplier);
@@ -1255,6 +1262,12 @@
 }
 template <>
 EIGEN_DEVICE_FUNC inline Packet8d pgather<double, Packet8d>(const double* from, Index stride) {
+  if (stride == 2) {
+    const Packet8d low = _mm512_loadu_pd(from);
+    const Packet8d high = _mm512_loadu_pd(from + 7);
+    const __m512i indices = _mm512_setr_epi64(0, 2, 4, 6, 9, 11, 13, 15);
+    return _mm512_permutex2var_pd(low, indices, high);
+  }
   Packet8i stride_vector = _mm256_set1_epi32(convert_index<int>(stride));
   Packet8i stride_multiplier = _mm256_set_epi32(7, 6, 5, 4, 3, 2, 1, 0);
   Packet8i indices = _mm256_mullo_epi32(stride_vector, stride_multiplier);

diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h
index 4e91e4d..ff3d13a 100644
--- a/Eigen/src/Core/arch/SSE/PacketMath.h
+++ b/Eigen/src/Core/arch/SSE/PacketMath.h

@@ -1567,6 +1567,10 @@
 
 template <>
 EIGEN_STRONG_INLINE Packet4f pgather<float, Packet4f>(const float* from, Index stride) {
+  if (stride == 2) {
+    // Overlap the loads so the second ends at the last gathered coefficient, from[6].
+    return _mm_shuffle_ps(_mm_loadu_ps(from), _mm_loadu_ps(from + 3), _MM_SHUFFLE(3, 1, 2, 0));
+  }
   return _mm_set_ps(from[3 * stride], from[2 * stride], from[1 * stride], from[0 * stride]);
 }
 template <>

diff --git a/benchmarks/Core/CMakeLists.txt b/benchmarks/Core/CMakeLists.txt
index 67db020..1da6867 100644
--- a/benchmarks/Core/CMakeLists.txt
+++ b/benchmarks/Core/CMakeLists.txt

@@ -22,6 +22,8 @@
 eigen_add_benchmark(bench_move_semantics bench_move_semantics.cpp)
 eigen_add_benchmark(bench_reductions bench_reductions.cpp)
 eigen_add_benchmark(bench_vectorwise_norm bench_vectorwise_norm.cpp)
+eigen_add_benchmark(bench_complex_components bench_complex_components.cpp)
+eigen_add_benchmark(bench_pgather bench_pgather.cpp)
 eigen_add_benchmark(bench_boolean_reductions bench_boolean_reductions.cpp)
 eigen_add_benchmark(bench_count bench_count.cpp)
 eigen_add_benchmark(bench_dot bench_dot.cpp)

diff --git a/benchmarks/Core/bench_complex_components.cpp b/benchmarks/Core/bench_complex_components.cpp
new file mode 100644
index 0000000..3b5c6b8
--- /dev/null
+++ b/benchmarks/Core/bench_complex_components.cpp

@@ -0,0 +1,51 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+#include <Eigen/Core>
+
+using namespace Eigen;
+
+template <bool Imaginary, int NaNPropagation, typename Vector>
+EIGEN_DONT_INLINE typename Vector::RealScalar component_abs_max(Vector& values) {
+  if (Imaginary) return values.imag().cwiseAbs().template maxCoeff<NaNPropagation>();
+  return values.real().cwiseAbs().template maxCoeff<NaNPropagation>();
+}
+
+template <typename Real, bool Imaginary, bool Const, int NaNPropagation>
+static void BM_ComponentAbsMax(benchmark::State& state) {
+  using Complex = std::complex<Real>;
+  using Vector = Matrix<Complex, Dynamic, 1>;
+  Vector values = Vector::Random(state.range(0));
+  std::conditional_t<Const, const Vector, Vector>& input = values;
+  Real expected = 0;
+  for (Index i = 0; i < values.size(); ++i) {
+    const Real value = numext::abs(Imaginary ? values(i).imag() : values(i).real());
+    if (value > expected) expected = value;
+  }
+  if (component_abs_max<Imaginary, NaNPropagation>(input) != expected) {
+    state.SkipWithError("Incorrect component maximum");
+    return;
+  }
+  for (auto _ : state) {
+    benchmark::ClobberMemory();
+    benchmark::DoNotOptimize(component_abs_max<Imaginary, NaNPropagation>(input));
+  }
+  state.SetBytesProcessed(state.iterations() * values.size() * sizeof(Complex));
+}
+
+// Small packet boundaries, the issue's cache-resident workload, and larger working sets.
+#define COMPONENT_SIZES ->Arg(3)->Arg(8)->Arg(9)->Arg(64)->Arg(4096)->Arg(65536)
+#define REGISTER_COMPONENT(Real, Imaginary, Const, Policy) \
+  BENCHMARK_TEMPLATE(BM_ComponentAbsMax, Real, Imaginary, Const, Policy) COMPONENT_SIZES
+#define REGISTER_TYPE(Real)                              \
+  REGISTER_COMPONENT(Real, false, false, PropagateFast); \
+  REGISTER_COMPONENT(Real, false, false, PropagateNaN);  \
+  REGISTER_COMPONENT(Real, false, true, PropagateNaN);   \
+  REGISTER_COMPONENT(Real, true, false, PropagateNaN);   \
+  REGISTER_COMPONENT(Real, true, true, PropagateNaN)
+REGISTER_TYPE(float);
+REGISTER_TYPE(double);
+#undef REGISTER_TYPE
+#undef REGISTER_COMPONENT
+#undef COMPONENT_SIZES

diff --git a/benchmarks/Core/bench_pgather.cpp b/benchmarks/Core/bench_pgather.cpp
new file mode 100644
index 0000000..38c7894
--- /dev/null
+++ b/benchmarks/Core/bench_pgather.cpp

@@ -0,0 +1,43 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+#include <Eigen/Core>
+
+using namespace Eigen;
+
+template <typename Scalar, int Stride>
+EIGEN_DONT_INLINE void gather_packets(const Scalar* input, Scalar* output, Index size, Index runtime_stride) {
+  using Packet = typename internal::packet_traits<Scalar>::type;
+  constexpr Index PacketSize = internal::packet_traits<Scalar>::size;
+  const Index stride = Stride == Dynamic ? runtime_stride : Stride;
+  for (Index i = 0; i < size; i += PacketSize) {
+    internal::pstoreu(output + i, internal::pgather<Scalar, Packet>(input + i * stride, stride));
+  }
+}
+
+template <typename Scalar, int Stride>
+static void BM_Pgather(benchmark::State& state) {
+  const Index size = state.range(0);
+  const Index stride = Stride == Dynamic ? state.range(1) : Stride;
+  Array<Scalar, Dynamic, 1> input((size - 1) * stride + 1), output(size);
+  for (Index i = 0; i < input.size(); ++i) input(i) = Scalar(i % 127);
+  gather_packets<Scalar, Stride>(input.data(), output.data(), size, stride);
+  for (Index i = 0; i < size; ++i) {
+    if (output(i) != input(i * stride)) {
+      state.SkipWithError("Incorrect gathered coefficient");
+      return;
+    }
+  }
+  for (auto _ : state) {
+    benchmark::ClobberMemory();
+    gather_packets<Scalar, Stride>(input.data(), output.data(), size, stride);
+    benchmark::DoNotOptimize(output.data());
+  }
+  state.SetItemsProcessed(state.iterations() * size);
+}
+
+BENCHMARK_TEMPLATE(BM_Pgather, float, 2)->Arg(64)->Arg(4096);
+BENCHMARK_TEMPLATE(BM_Pgather, double, 2)->Arg(64)->Arg(4096);
+BENCHMARK_TEMPLATE(BM_Pgather, float, Dynamic)->ArgsProduct({{64, 4096}, {1, 2, 3, 7}});
+BENCHMARK_TEMPLATE(BM_Pgather, double, Dynamic)->ArgsProduct({{64, 4096}, {1, 2, 3, 7}});

diff --git a/benchmarks/Core/bench_reductions.cpp b/benchmarks/Core/bench_reductions.cpp
index e09f69e..16459d6 100644
--- a/benchmarks/Core/bench_reductions.cpp
+++ b/benchmarks/Core/bench_reductions.cpp

@@ -79,8 +79,7 @@
   state.SetBytesProcessed(state.iterations() * n * sizeof(Scalar));
 }
 
-// real() builds a CwiseUnaryView, which drops PacketAccessBit, so this reduces coefficient by
-// coefficient and reaches the scalar form of the wrapper.
+// Component reductions can use stride-2 packet gathers.
 template <typename Scalar, int NaNPropagation>
 static void BM_ComplexRealAbsMaxCoeff(benchmark::State& state) {
   using Real = typename NumTraits<Scalar>::Real;

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 79da7f4..553fb80 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt

@@ -269,6 +269,7 @@
 ei_add_test(clz)
 ei_add_test(rand)
 ei_add_test(realview)
+ei_add_test(complex_components)
 ei_add_test(meta)
 ei_add_test(preprocessor_directives)
 ei_add_test(lru_cache)
@@ -292,6 +293,7 @@
 ei_add_test(float_conversion)
 ei_add_test(io)
 ei_add_test(packetmath "-DEIGEN_FAST_MATH=1")
+ei_add_test(packetmath_gather)
 # Generic clang vector backend tests for different vector sizes.
 include(CheckCXXSourceCompiles)
 check_cxx_source_compiles("

diff --git a/test/complex_components.cpp b/test/complex_components.cpp
new file mode 100644
index 0000000..c6bef59
--- /dev/null
+++ b/test/complex_components.cpp

@@ -0,0 +1,177 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include "main.h"
+#include <Eigen/Core>
+
+// Component access must not assume that a user-defined complex scalar has std::complex's layout.
+struct PaddedComplex : std::complex<double> {
+  using std::complex<double>::complex;
+  double padding = 123;
+};
+
+template <>
+struct Eigen::NumTraits<PaddedComplex> : Eigen::NumTraits<std::complex<double>> {};
+
+template <typename Derived>
+void check_component(const DenseBase<Derived>& component) {
+  using Real = typename Derived::Scalar;
+  Matrix<Real, Dynamic, Dynamic, Derived::IsRowMajor ? RowMajor : ColMajor> copy = component;
+  Real sum = 0, minimum = component(0, 0), maximum = minimum, abs_maximum = numext::abs(minimum);
+  for (Index col = 0; col < component.cols(); ++col) {
+    for (Index row = 0; row < component.rows(); ++row) {
+      const Real value = component(row, col);
+      VERIFY_IS_EQUAL(copy(row, col), value);
+      sum += value;
+      if (value < minimum) minimum = value;
+      if (value > maximum) maximum = value;
+      if (numext::abs(value) > abs_maximum) abs_maximum = numext::abs(value);
+    }
+  }
+  VERIFY_IS_EQUAL(component.sum(), sum);
+  VERIFY_IS_EQUAL(component.template minCoeff<PropagateNaN>(), minimum);
+  VERIFY_IS_EQUAL(component.template maxCoeff<PropagateNaN>(), maximum);
+  VERIFY_IS_EQUAL(component.derived().cwiseAbs().template maxCoeff<PropagateNaN>(), abs_maximum);
+  VERIFY_IS_CWISE_EQUAL(component.derived().colwise().sum().eval(), copy.colwise().sum().eval());
+  VERIFY_IS_CWISE_EQUAL(component.derived().rowwise().sum().eval(), copy.rowwise().sum().eval());
+}
+
+template <typename MatrixType>
+void component_expressions(Index rows, Index cols) {
+  using Complex = typename MatrixType::Scalar;
+  using Real = typename NumTraits<Complex>::Real;
+  MatrixType values(rows, cols);
+  for (Index i = 0; i < values.size(); ++i) values.data()[i] = Complex(Real(i % 13 - 6), Real(i % 17 + 1));
+  const MatrixType& const_values = values;
+  check_component(values.real());
+  check_component(values.imag());
+  check_component(const_values.real());
+  check_component(const_values.imag());
+  check_component(values.transpose().real());
+  check_component(values.transpose().imag());
+  check_component(values.conjugate().real());
+  check_component(values.conjugate().imag());
+  check_component(values.block(0, 0, rows - 1, cols - 1).real());
+  check_component(values.block(0, 0, rows - 1, cols - 1).imag());
+  check_component(const_values.block(0, 0, rows - 1, cols - 1).real());
+  check_component(const_values.block(0, 0, rows - 1, cols - 1).imag());
+  const auto real = values.real();
+  const auto imag = values.imag();
+  check_component(real.block(0, 0, rows - 1, cols - 1));
+  check_component(imag.block(0, 0, rows - 1, cols - 1));
+
+  MatrixType original = values;
+  values.real() = original.imag();
+  values.imag() = original.real();
+  for (Index i = 0; i < values.size(); ++i) {
+    VERIFY_IS_EQUAL(values.data()[i], Complex(original.data()[i].imag(), original.data()[i].real()));
+  }
+  values.real().block(0, 0, rows - 1, cols - 1).setConstant(Real(3));
+  values.imag().array() += Real(2);
+  for (Index col = 0; col < cols; ++col) {
+    for (Index row = 0; row < rows; ++row) {
+      VERIFY_IS_EQUAL(values(row, col).imag(), original(row, col).real() + Real(2));
+      VERIFY_IS_EQUAL(values(row, col).real(), row < rows - 1 && col < cols - 1 ? Real(3) : original(row, col).imag());
+    }
+  }
+}
+
+template <typename Real>
+void component_reductions() {
+  using Complex = std::complex<Real>;
+  using Vector = Matrix<Complex, Dynamic, 1>;
+  constexpr bool vectorizable = (std::is_same<Real, float>::value || std::is_same<Real, double>::value) &&
+                                internal::packet_traits<Real>::Vectorizable;
+  using View = decltype(std::declval<Vector&>().real());
+  using ConstOp = decltype(std::declval<const Vector&>().imag());
+  STATIC_CHECK(bool(internal::redux_evaluator<View>::Flags & PacketAccessBit) == vectorizable);
+  STATIC_CHECK(bool(internal::evaluator<ConstOp>::Flags & PacketAccessBit) == vectorizable);
+  STATIC_CHECK(!(internal::evaluator<View>::Flags & PacketAccessBit));
+  using Strided = Map<Vector, Unaligned, InnerStride<2>>;
+  using StridedView = decltype(std::declval<Strided&>().real());
+  STATIC_CHECK(!(internal::redux_evaluator<StridedView>::Flags & PacketAccessBit));
+  using DynamicStride = Map<Vector, Unaligned, InnerStride<Dynamic>>;
+  using DynamicView = decltype(std::declval<const DynamicStride&>().imag());
+  STATIC_CHECK(!(internal::evaluator<DynamicView>::Flags & PacketAccessBit));
+  using SumView = decltype((std::declval<Vector>() + std::declval<Vector>()).real());
+  STATIC_CHECK(!(internal::evaluator<SumView>::Flags & PacketAccessBit));
+
+  const Real nan = NumTraits<Real>::quiet_NaN();
+  const Real inf = NumTraits<Real>::infinity();
+  constexpr Index packet_size = internal::packet_traits<Real>::size;
+  for (Index size : {Index(1), packet_size, packet_size + 1, 2 * packet_size - 1, Index(65)}) {
+    Vector buffer(2 * size + 1);
+    buffer.setConstant(Complex(Real(7), Real(-3)));
+    Map<Vector> values(buffer.data() + 1, size);
+    const auto& const_values = values;
+    check_component(values.real());
+    check_component(values.imag());
+    Strided strided(buffer.data() + 1, size);
+    check_component(strided.real());
+    check_component(strided.imag());
+    DynamicStride dynamic(buffer.data() + 1, size, InnerStride<Dynamic>(2));
+    check_component(dynamic.real());
+    check_component(dynamic.imag());
+    for (Index i = 0; i < size; ++i) {
+      for (Real special : {nan, inf, -inf, Real(0), -Real(0), std::numeric_limits<Real>::denorm_min()}) {
+        values.setConstant(Complex(Real(7), Real(-3)));
+        values(i) = Complex(special, Real(-3));
+        Matrix<Real, Dynamic, 1> real = const_values.real();
+        if ((numext::isnan)(special)) {
+          VERIFY((numext::isnan)(real(i)));
+          VERIFY((numext::isnan)(values.real().template maxCoeff<PropagateNaN>()));
+          VERIFY((numext::isnan)(const_values.real().template minCoeff<PropagateNaN>()));
+          if (size > 1) VERIFY_IS_EQUAL(values.real().template maxCoeff<PropagateNumbers>(), Real(7));
+        } else {
+          VERIFY_IS_EQUAL(real(i), special);
+          VERIFY_IS_EQUAL((std::signbit)(real(i)), (std::signbit)(special));
+        }
+        VERIFY_IS_EQUAL(values.imag().template maxCoeff<PropagateNaN>(), Real(-3));
+        values(i) = Complex(Real(7), special);
+        Matrix<Real, Dynamic, 1> imag = const_values.imag();
+        if ((numext::isnan)(special)) {
+          VERIFY((numext::isnan)(imag(i)));
+          VERIFY((numext::isnan)(values.imag().template minCoeff<PropagateNaN>()));
+          VERIFY((numext::isnan)(const_values.imag().template maxCoeff<PropagateNaN>()));
+          if (size > 1) VERIFY_IS_EQUAL(values.imag().template minCoeff<PropagateNumbers>(), Real(-3));
+        } else {
+          VERIFY_IS_EQUAL(imag(i), special);
+          VERIFY_IS_EQUAL((std::signbit)(imag(i)), (std::signbit)(special));
+        }
+        VERIFY_IS_EQUAL(values.real().template minCoeff<PropagateNaN>(), Real(7));
+      }
+    }
+  }
+}
+
+void component_custom_scalar() {
+  using Vector = Matrix<PaddedComplex, Dynamic, 1>;
+  Vector values(17);
+  for (Index i = 0; i < values.size(); ++i) values(i) = PaddedComplex(double(i), double(2 * i + 1));
+  const Vector& const_values = values;
+  using RealView = decltype(const_values.real());
+  using ImagView = decltype(const_values.imag());
+  STATIC_CHECK(!(internal::evaluator<RealView>::Flags & PacketAccessBit));
+  STATIC_CHECK(!(internal::evaluator<ImagView>::Flags & PacketAccessBit));
+  check_component(const_values.real());
+  check_component(const_values.imag());
+}
+
+EIGEN_DECLARE_TEST(complex_components) {
+  CALL_SUBTEST_1(component_reductions<float>());
+  CALL_SUBTEST_2(component_reductions<double>());
+  CALL_SUBTEST_3(component_reductions<long double>());
+  CALL_SUBTEST_3(component_custom_scalar());
+  for (Index rows : {Index(2), Index(7), Index(17)}) {
+    for (Index cols : {Index(3), Index(9)}) {
+      EIGEN_UNUSED_VARIABLE(rows);
+      EIGEN_UNUSED_VARIABLE(cols);
+      CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, Dynamic, Dynamic, ColMajor>>(rows, cols)));
+      CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>>(rows, cols)));
+      CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, Dynamic, Dynamic, ColMajor>>(rows, cols)));
+      CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>(rows, cols)));
+    }
+  }
+  CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, 3, 5>>(3, 5)));
+  CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, 2, 2>>(2, 2)));
+}

diff --git a/test/packetmath_gather.cpp b/test/packetmath_gather.cpp
new file mode 100644
index 0000000..f88e6c4
--- /dev/null
+++ b/test/packetmath_gather.cpp

@@ -0,0 +1,48 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include "packetmath_test_shared.h"
+#include <memory>
+
+namespace Eigen {
+namespace test {
+
+template <typename Scalar, typename Packet>
+struct runall<Scalar, Packet, false, false> {
+  static void run() {
+    constexpr int PacketSize = internal::unpacket_traits<Packet>::size;
+    auto samples = special_values<Scalar>();
+    using Bits = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type;
+    const Bits nan_payload = numext::bit_cast<Bits>(NumTraits<Scalar>::quiet_NaN()) | Bits(0x123);
+    samples.push_back(numext::bit_cast<Scalar>(nan_payload));
+    Scalar actual[PacketSize], expected[PacketSize];
+    for (Index stride : {Index(-7), Index(-2), Index(-1), Index(0), Index(1), Index(2), Index(3), Index(7)}) {
+      const Index span = (PacketSize - 1) * numext::abs(stride) + 1;
+      for (Index offset = 0; offset < PacketSize; ++offset) {
+        // Exact-size new[] lets sanitizers detect reads past the last gathered coefficient.
+        std::unique_ptr<Scalar[]> storage(new Scalar[span + offset]);
+        Scalar* data = storage.get() + offset;
+        const Scalar* from = stride < 0 ? data + span - 1 : data;
+        for (Index rotation = 0; rotation < Index(samples.size()); ++rotation) {
+          for (Index i = 0; i < span; ++i) data[i] = samples[(i + rotation) % samples.size()];
+          for (Index i = 0; i < PacketSize; ++i) expected[i] = from[i * stride];
+          volatile Index runtime_stride = stride;
+          internal::pstoreu(actual, internal::pgather<Scalar, Packet>(from, runtime_stride));
+          VERIFY(areEqualBits(expected, actual, PacketSize, false));
+          if (stride == 2) {
+            internal::pstoreu(actual, internal::pgather<Scalar, Packet>(from, 2));
+            VERIFY(areEqualBits(expected, actual, PacketSize, false));
+          }
+        }
+      }
+    }
+  }
+};
+
+}  // namespace test
+}  // namespace Eigen
+
+EIGEN_DECLARE_TEST(packetmath_gather) {
+  CALL_SUBTEST_1(test::runner<float>::run());
+  CALL_SUBTEST_2(test::runner<double>::run());
+}