Core: Work around GCC PR124138 in unrolled self-adjoint copies libeigen/eigen!3411 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
index 21cffb3..b1bcf57 100644
--- a/Eigen/src/Core/SelfAdjointView.h
+++ b/Eigen/src/Core/SelfAdjointView.h

@@ -379,6 +379,23 @@
   using Shape = SelfAdjointShape;
 };
 
+// GCC PR target/124138 (https://gcc.gnu.org/PR124138): since GCC 4.9, avx_vpermilp_parallel accepted selectors 4-7
+// for elements 2 and 3 of a V8DF permutation, emitting the lane-crossing {0, 1, 2, 6, 4, 5, 6, 7} as vpermilpd, which
+// yields {0, 1, 2, 2, 4, 5, 6, 7}. SLP forms that permutation from the unrolled mirrored stores of a 4x4 double copy,
+// at -O2 since GCC 12. Hiding the copied value from the vectorizer avoids it; the loop used for larger sizes is
+// unaffected. Fixed in GCC 14.4, 15.3 and 16, and on the GCC 13 branch after 13.4.
+template <bool Unrolled, typename Scalar>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void selfadjoint_mirror_barrier(Scalar&) {}
+#if EIGEN_COMP_GNUC_STRICT && defined(__AVX512F__) && !defined(EIGEN_GPU_COMPILE_PHASE) && \
+    (EIGEN_GNUC_STRICT_LESS_THAN(13, 5, 0) ||                                              \
+     (EIGEN_GNUC_STRICT_AT_LEAST(14, 0, 0) && EIGEN_GNUC_STRICT_LESS_THAN(14, 4, 0)) ||    \
+     (EIGEN_GNUC_STRICT_AT_LEAST(15, 0, 0) && EIGEN_GNUC_STRICT_LESS_THAN(15, 3, 0)))
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void selfadjoint_mirror_barrier<true, double>(double& value) {
+  EIGEN_OPTIMIZATION_BARRIER(value);
+}
+#endif
+
 template <int UpLo, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor,
           int Version>
 class triangular_dense_assignment_kernel<UpLo, SelfAdjoint, SetOpposite, DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor,
@@ -405,6 +422,7 @@
   EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col) {
     eigen_internal_assert(row != col);
     Scalar tmp = m_src.coeff(row, col);
+    selfadjoint_mirror_barrier<triangular_assignment_unrolls<DstEvaluatorType, SrcEvaluatorType>::value>(tmp);
     m_functor.assignCoeff(m_dst.coeffRef(row, col), tmp);
     m_functor.assignCoeff(m_dst.coeffRef(col, row), numext::conj(tmp));
   }

diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index 0d083aa..ee0e14f 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h

@@ -846,6 +846,15 @@
 template <typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
 struct triangular_assignment_loop;
 
+// Whether call_triangular_assignment_loop fully unrolls the assignment.
+template <typename DstEvaluatorType, typename SrcEvaluatorType>
+struct triangular_assignment_unrolls {
+  static constexpr int Size = DstEvaluatorType::XprType::SizeAtCompileTime;
+  static constexpr bool value =
+      Size != Dynamic && SrcEvaluatorType::CoeffReadCost < HugeCost &&
+      Size * (int(DstEvaluatorType::CoeffReadCost) + int(SrcEvaluatorType::CoeffReadCost)) / 2 <= EIGEN_UNROLLING_LIMIT;
+};
+
 /** \internal Specialization of the dense assignment kernel for triangular matrices.
  * The main difference is that the triangular, diagonal, and opposite parts are processed through three different
  * functions. \tparam UpLo must be either Lower or Upper \tparam Mode must be either 0, UnitDiag, ZeroDiag, or
@@ -917,13 +926,7 @@
                                                     SetOpposite, DstEvaluatorType, SrcEvaluatorType, Functor>;
   Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
 
-  enum {
-    unroll = DstXprType::SizeAtCompileTime != Dynamic && SrcEvaluatorType::CoeffReadCost < HugeCost &&
-             DstXprType::SizeAtCompileTime *
-                     (int(DstEvaluatorType::CoeffReadCost) + int(SrcEvaluatorType::CoeffReadCost)) / 2 <=
-                 EIGEN_UNROLLING_LIMIT
-  };
-
+  constexpr bool unroll = triangular_assignment_unrolls<DstEvaluatorType, SrcEvaluatorType>::value;
   triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(
       kernel);
 }

diff --git a/test/selfadjoint.cpp b/test/selfadjoint.cpp
index 359a8a0..109a719 100644
--- a/test/selfadjoint.cpp
+++ b/test/selfadjoint.cpp

@@ -175,6 +175,28 @@
   EIGEN_UNUSED_VARIABLE(m);
 }
 
+template <unsigned int UpLo, typename MatrixType>
+EIGEN_DONT_INLINE void verify_selfadjoint_copy(const MatrixType& m, const MatrixType& dense) {
+  for (Index j = 0; j < m.cols(); ++j)
+    for (Index i = 0; i < m.rows(); ++i)
+      VERIFY_IS_EQUAL(dense(i, j), (UpLo == Upper ? i <= j : i >= j) ? m(i, j) : numext::conj(m(j, i)));
+}
+
+// GCC PR target/124138 miscompiled this copy for 4x4 double with AVX-512 (see SelfAdjointView.h). It needs the copy
+// materialized in a function of its own, and a reference built without the mirrored stores that GCC also miscompiles.
+template <unsigned int UpLo, typename MatrixType>
+EIGEN_DONT_INLINE void selfadjoint_copy(const MatrixType& m) {
+  const MatrixType dense = m.template selfadjointView<UpLo>();
+  verify_selfadjoint_copy<UpLo>(m, dense);
+}
+
+template <typename MatrixType>
+void selfadjoint_copy_fixed() {
+  const MatrixType m = MatrixType::Random();
+  selfadjoint_copy<Upper>(m);
+  selfadjoint_copy<Lower>(m);
+}
+
 EIGEN_DECLARE_TEST(selfadjoint) {
   for (int i = 0; i < g_repeat; i++) {
     int s = internal::random<int>(1, EIGEN_TEST_MAX_SIZE);
@@ -190,6 +212,9 @@
   }
 
   CALL_SUBTEST_1(bug_159());
+  CALL_SUBTEST_2(selfadjoint_copy_fixed<Matrix4d>());
+  CALL_SUBTEST_2((selfadjoint_copy_fixed<Matrix<double, 4, 4, RowMajor> >()));
+  CALL_SUBTEST_3(selfadjoint_copy_fixed<Matrix4cd>());
   CALL_SUBTEST_1(selfadjoint_nested_expression<>());
   CALL_SUBTEST_4(selfadjoint_l1norm_sizes<double>());
   CALL_SUBTEST_4(selfadjoint_l1norm_sizes<std::complex<double> >());