SelfAdjointEigenSolver: scalar argmin in eigenvalue sort (fixes ARM NEON sort)

libeigen/eigen!2505

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index 65533d8..c89c2d1 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -626,11 +626,21 @@
   // TODO: make the sort optional and use a more efficient sorting algorithm.
   if (info == Success) {
     for (Index i = 0; i < n - 1; ++i) {
-      Index k;
-      diag.segment(i, n - i).minCoeff(&k);
-      if (k > 0) {
-        numext::swap(diag[i], diag[k + i]);
-        if (computeEigenvectors) eivec.col(i).swap(eivec.col(k + i));
+      // Scalar argmin: the vectorized minCoeff path can return the wrong
+      // index on targets whose SIMD reduction flushes subnormals (32-bit ARM
+      // NEON always treats subnormal inputs to vminq_f32 as zero), which
+      // corrupts ordering when eigenvalues span the subnormal range.
+      Index k = i;
+      RealScalar min_val = diag[i];
+      for (Index j = i + 1; j < n; ++j) {
+        if (diag[j] < min_val) {
+          min_val = diag[j];
+          k = j;
+        }
+      }
+      if (k != i) {
+        numext::swap(diag[i], diag[k]);
+        if (computeEigenvectors) eivec.col(i).swap(eivec.col(k));
       }
     }
   }