From be20207d108a9d74d62bfe3e7d8680a17483551e Mon Sep 17 00:00:00 2001
From: Arthur <afeeney@umass.edu>
Date: Mon, 8 Aug 2022 19:29:56 +0000
Subject: [PATCH] Fix vectorized Jacobi Rotation

---
 Eigen/src/Jacobi/Jacobi.h | 23 +++++++++++++----------
 test/jacobi.cpp           |  5 +++++
 2 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h
index 852a38533..5d9698990 100644
--- a/Eigen/src/Jacobi/Jacobi.h
+++ b/Eigen/src/Jacobi/Jacobi.h
@@ -344,15 +344,18 @@ struct apply_rotation_in_the_plane_selector<Scalar,OtherScalar,SizeAtCompileTime
 {
   static inline void run(Scalar *x, Index incrx, Scalar *y, Index incry, Index size, OtherScalar c, OtherScalar s)
   {
-    enum {
-      PacketSize = packet_traits<Scalar>::size,
-      OtherPacketSize = packet_traits<OtherScalar>::size
-    };
     typedef typename packet_traits<Scalar>::type Packet;
     typedef typename packet_traits<OtherScalar>::type OtherPacket;
 
+    enum {
+      RequiredAlignment = plain_enum_max(unpacket_traits<Packet>::alignment,
+                                         unpacket_traits<OtherPacket>::alignment),
+      PacketSize = packet_traits<Scalar>::size,
+      OtherPacketSize = packet_traits<OtherScalar>::size
+    };
+
     /*** dynamic-size vectorized paths ***/
-    if(SizeAtCompileTime == Dynamic && ((incrx==1 && incry==1) || PacketSize == 1))
+    if(size >= 2 * PacketSize && SizeAtCompileTime == Dynamic && ((incrx == 1 && incry == 1) || PacketSize == 1))
     {
       // both vectors are sequentially stored in memory => vectorization
       enum { Peeling = 2 };
@@ -423,11 +426,11 @@ struct apply_rotation_in_the_plane_selector<Scalar,OtherScalar,SizeAtCompileTime
     }
 
     /*** fixed-size vectorized path ***/
-    else if(SizeAtCompileTime != Dynamic && MinAlignment>0) // FIXME should be compared to the required alignment
+    else if(SizeAtCompileTime != Dynamic && MinAlignment >= RequiredAlignment)
     {
       const OtherPacket pc = pset1<OtherPacket>(c);
       const OtherPacket ps = pset1<OtherPacket>(s);
-      conj_helper<OtherPacket,Packet,NumTraits<OtherPacket>::IsComplex,false> pcj;
+      conj_helper<OtherPacket,Packet,NumTraits<OtherScalar>::IsComplex,false> pcj;
       conj_helper<OtherPacket,Packet,false,false> pm;
       Scalar* EIGEN_RESTRICT px = x;
       Scalar* EIGEN_RESTRICT py = y;
@@ -452,11 +455,11 @@ struct apply_rotation_in_the_plane_selector<Scalar,OtherScalar,SizeAtCompileTime
 
 template<typename VectorX, typename VectorY, typename OtherScalar>
 EIGEN_DEVICE_FUNC
-void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(DenseBase<VectorX>& xpr_x, DenseBase<VectorY>& xpr_y, const JacobiRotation<OtherScalar>& j)
+void inline apply_rotation_in_the_plane(DenseBase<VectorX>& xpr_x, DenseBase<VectorY>& xpr_y, const JacobiRotation<OtherScalar>& j)
 {
   typedef typename VectorX::Scalar Scalar;
-  const bool Vectorizable =    (int(VectorX::Flags) & int(VectorY::Flags) & PacketAccessBit)
-                            && (int(packet_traits<Scalar>::size) == int(packet_traits<OtherScalar>::size));
+  constexpr bool Vectorizable = (int(evaluator<VectorX>::Flags) & int(evaluator<VectorY>::Flags) & PacketAccessBit) &&
+                                (int(packet_traits<Scalar>::size) == int(packet_traits<OtherScalar>::size));
 
   eigen_assert(xpr_x.size() == xpr_y.size());
   Index size = xpr_x.size();
diff --git a/test/jacobi.cpp b/test/jacobi.cpp
index 5604797f5..273b94df8 100644
--- a/test/jacobi.cpp
+++ b/test/jacobi.cpp
@@ -65,6 +65,11 @@ EIGEN_DECLARE_TEST(jacobi)
     CALL_SUBTEST_3(( jacobi<Matrix4cf, float>() ));
     CALL_SUBTEST_3(( jacobi<Matrix4cf, std::complex<float> >() ));
 
+    CALL_SUBTEST_1(( jacobi<Matrix<float, 3, 3, RowMajor>, float>() ));
+    CALL_SUBTEST_2(( jacobi<Matrix<double, 4, 4, RowMajor>, double>() ));
+    CALL_SUBTEST_3(( jacobi<Matrix<std::complex<float>, 4, 4, RowMajor>, float>() ));
+    CALL_SUBTEST_3(( jacobi<Matrix<std::complex<float>, 4, 4, RowMajor>, std::complex<float> >() ));
+
     int r = internal::random<int>(2, internal::random<int>(1,EIGEN_TEST_MAX_SIZE)/2),
         c = internal::random<int>(2, internal::random<int>(1,EIGEN_TEST_MAX_SIZE)/2);
     CALL_SUBTEST_4(( jacobi<MatrixXf, float>(MatrixXf(r,c)) ));