From ea13a98decd497a8c5588fb5de71b57bcf10d864 Mon Sep 17 00:00:00 2001
From: Rasmus Munk Larsen <4643818-rmlarsen1@users.noreply.gitlab.com>
Date: Sun, 15 Mar 2026 19:56:01 -0700
Subject: [PATCH] Fix imag_ref for real scalar types and clean up svd_fill.h

libeigen/eigen!2303

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
---
 Eigen/src/Core/MathFunctions.h | 12 ++++++---
 test/svd_fill.h                | 47 ++++++++++++++++++++++++++--------
 2 files changed, 45 insertions(+), 14 deletions(-)

diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index 4e7597ced..54da17cc1 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -170,18 +170,24 @@ struct imag_ref_default_impl {
 
 template <typename Scalar>
 struct imag_ref_default_impl<Scalar, false> {
-  EIGEN_DEVICE_FUNC constexpr static Scalar run(Scalar&) { return Scalar(0); }
-  EIGEN_DEVICE_FUNC constexpr static const Scalar run(const Scalar&) { return Scalar(0); }
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  EIGEN_DEVICE_FUNC constexpr static inline RealScalar run(Scalar&) { return RealScalar(0); }
+  EIGEN_DEVICE_FUNC constexpr static inline RealScalar run(const Scalar&) { return RealScalar(0); }
 };
 
 template <typename Scalar>
 struct imag_ref_impl : imag_ref_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
 
-template <typename Scalar>
+template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
 struct imag_ref_retval {
   typedef typename NumTraits<Scalar>::Real& type;
 };
 
+template <typename Scalar>
+struct imag_ref_retval<Scalar, false> {
+  typedef typename NumTraits<Scalar>::Real type;
+};
+
 }  // namespace internal
 
 namespace numext {
diff --git a/test/svd_fill.h b/test/svd_fill.h
index d092e836b..3ef7c9ae3 100644
--- a/test/svd_fill.h
+++ b/test/svd_fill.h
@@ -23,6 +23,15 @@ Array<T, 4, 1> four_denorms() {
   return four_denorms<double>().cast<T>();
 }
 
+template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
+struct maybe_set_imag_part {
+  static void run(Scalar& x, typename NumTraits<Scalar>::Real y) { numext::imag_ref(x) = y; }
+};
+template <typename Scalar>
+struct maybe_set_imag_part<Scalar, false> {
+  static void run(Scalar&, typename NumTraits<Scalar>::Real) {}
+};
+
 template <typename MatrixType>
 void svd_fill_random(MatrixType &m, int Option = 0) {
   using std::pow;
@@ -35,9 +44,12 @@ void svd_fill_random(MatrixType &m, int Option = 0) {
   for (Index k = 0; k < diagSize; ++k) d(k) = d(k) * pow(RealScalar(10), internal::random<RealScalar>(-s, s));
 
   bool dup = internal::random<int>(0, 10) < 3;
-  bool unit_uv =
-      internal::random<int>(0, 10) < (dup ? 7 : 3);  // if we duplicate some diagonal entries, then increase the chance
-                                                     // to preserve them using unitary U and V factors
+  // For SelfAdjoint, always use unitary U so the eigenvalues are exactly d
+  // and the conditioning is controlled.  For Symmetric/general SVD, randomly
+  // choose unitary or arbitrary U/V.
+  bool unit_uv = (Option == SelfAdjoint) ||
+                 internal::random<int>(0, 10) < (dup ? 7 : 3);  // if we duplicate some diagonal entries, then increase
+                                                                // the chance to preserve them using unitary U and V
 
   // duplicate some singular values
   if (dup) {
@@ -67,8 +79,11 @@ void svd_fill_random(MatrixType &m, int Option = 0) {
       RealScalar(1) / NumTraits<RealScalar>::highest(), (std::numeric_limits<RealScalar>::min)(),
       pow((std::numeric_limits<RealScalar>::min)(), RealScalar(0.8));
 
-  if (Option == Symmetric) {
-    m = U * d.asDiagonal() * U.transpose();
+  if (Option == Symmetric || Option == SelfAdjoint) {
+    if (Option == SelfAdjoint)
+      m = U * d.asDiagonal() * U.adjoint();
+    else
+      m = U * d.asDiagonal() * U.transpose();
 
     // randomly nullify some rows/columns
     {
@@ -85,10 +100,21 @@ void svd_fill_random(MatrixType &m, int Option = 0) {
           for (Index k = 0; k < n; ++k) {
             Index i = internal::random<Index>(0, m.rows() - 1);
             Index j = internal::random<Index>(0, m.cols() - 1);
-            m(j, i) = m(i, j) = samples(internal::random<Index>(0, samples.size() - 1));
-            if (NumTraits<Scalar>::IsComplex)
-              *(&numext::real_ref(m(j, i)) + 1) = *(&numext::real_ref(m(i, j)) + 1) =
-                  samples.real()(internal::random<Index>(0, samples.size() - 1));
+            Scalar val = samples(internal::random<Index>(0, samples.size() - 1));
+            maybe_set_imag_part<Scalar>::run(val, samples.real()(internal::random<Index>(0, samples.size() - 1)));
+            if (Option == SelfAdjoint) {
+              if (i == j) {
+                // Diagonal entries of a Hermitian matrix must be real.
+                m(i, i) = numext::real(val);
+              } else {
+                // Off-diagonal: m(j,i) = conj(m(i,j)) for Hermitian.
+                m(i, j) = val;
+                m(j, i) = numext::conj(val);
+              }
+            } else {
+              m(i, j) = val;
+              m(j, i) = val;
+            }
           }
         }
     }
@@ -101,8 +127,7 @@ void svd_fill_random(MatrixType &m, int Option = 0) {
         Index i = internal::random<Index>(0, m.rows() - 1);
         Index j = internal::random<Index>(0, m.cols() - 1);
         m(i, j) = samples(internal::random<Index>(0, samples.size() - 1));
-        if (NumTraits<Scalar>::IsComplex)
-          *(&numext::real_ref(m(i, j)) + 1) = samples.real()(internal::random<Index>(0, samples.size() - 1));
+        maybe_set_imag_part<Scalar>::run(m(i, j), samples.real()(internal::random<Index>(0, samples.size() - 1)));
       }
     }
   }