From cb8e6d497561c53523315e8d1931254a483dff5b Mon Sep 17 00:00:00 2001
From: Charles Schlosser <charless@vt.edu>
Date: Mon, 6 Mar 2023 23:11:06 +0000
Subject: [PATCH] Fix 2240, 2620

---
 .../src/Eigenvalues/SelfAdjointEigenSolver.h  |  7 ++++-
 Eigen/src/Eigenvalues/Tridiagonalization.h    | 28 ++++++++++---------
 2 files changed, 21 insertions(+), 14 deletions(-)

diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index d196ec054..88e05bd73 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -108,6 +108,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
       * This is a column vector with entries of type #RealScalar.
       * The length of the vector is the size of \p MatrixType_.
       */
+    typedef typename internal::plain_col_type<MatrixType, Scalar>::type VectorType;
     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
     typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType;
@@ -125,6 +126,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
     EIGEN_DEVICE_FUNC
     SelfAdjointEigenSolver()
         : m_eivec(),
+          m_workspace(),
           m_eivalues(),
           m_subdiag(),
           m_hcoeffs(),
@@ -148,6 +150,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
     EIGEN_DEVICE_FUNC
     explicit SelfAdjointEigenSolver(Index size)
         : m_eivec(size, size),
+          m_workspace(size),
           m_eivalues(size),
           m_subdiag(size > 1 ? size - 1 : 1),
           m_hcoeffs(size > 1 ? size - 1 : 1),
@@ -174,6 +177,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
     EIGEN_DEVICE_FUNC
     explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors)
       : m_eivec(matrix.rows(), matrix.cols()),
+        m_workspace(matrix.cols()),
         m_eivalues(matrix.cols()),
         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
         m_hcoeffs(matrix.cols() > 1 ? matrix.cols() - 1 : 1),
@@ -377,6 +381,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
     EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
 
     EigenvectorsType m_eivec;
+    VectorType m_workspace;
     RealVectorType m_eivalues;
     typename TridiagonalizationType::SubDiagonalType m_subdiag;
     typename TridiagonalizationType::CoeffVectorType m_hcoeffs;
@@ -451,7 +456,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
   mat.template triangularView<Lower>() /= scale;
   m_subdiag.resize(n-1);
   m_hcoeffs.resize(n-1);
-  internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, computeEigenvectors);
+  internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, m_workspace, computeEigenvectors);
 
   m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
   
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h
index 9b002fe7b..a975e563a 100644
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h
@@ -427,13 +427,13 @@ struct tridiagonalization_inplace_selector;
   *
   * \sa class Tridiagonalization
   */
-template<typename MatrixType, typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
+template<typename MatrixType, typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType, typename WorkSpaceType>
 EIGEN_DEVICE_FUNC
 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
-                                CoeffVectorType& hcoeffs, bool extractQ)
+                                CoeffVectorType& hcoeffs, WorkSpaceType& workspace, bool extractQ)
 {
   eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
-  tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, hcoeffs, extractQ);
+  tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, hcoeffs, workspace, extractQ);
 }
 
 /** \internal
@@ -443,17 +443,19 @@ template<typename MatrixType, int Size, bool IsComplex>
 struct tridiagonalization_inplace_selector
 {
   typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
-  template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
+  template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType, typename WorkSpaceType>
   static EIGEN_DEVICE_FUNC
-      void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType& hCoeffs, bool extractQ)
+      void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType& hCoeffs, WorkSpaceType& workspace, bool extractQ)
   {
     tridiagonalization_inplace(mat, hCoeffs);
     diag = mat.diagonal().real();
     subdiag = mat.template diagonal<-1>().real();
-    if(extractQ)
-      mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
-            .setLength(mat.rows() - 1)
-            .setShift(1);
+    if (extractQ) {
+      HouseholderSequenceType hh(mat, hCoeffs.conjugate());
+      hh.setLength(mat.rows() - 1);
+      hh.setShift(1);
+      hh.evalTo(mat, workspace);
+    }
   }
 };
 
@@ -467,8 +469,8 @@ struct tridiagonalization_inplace_selector<MatrixType,3,false>
   typedef typename MatrixType::Scalar Scalar;
   typedef typename MatrixType::RealScalar RealScalar;
 
-  template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
-  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType&, bool extractQ)
+  template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType, typename WorkSpaceType>
+  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType&, WorkSpaceType&, bool extractQ)
   {
     using std::sqrt;
     const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
@@ -512,9 +514,9 @@ struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
 {
   typedef typename MatrixType::Scalar Scalar;
 
-  template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
+  template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType, typename WorkSpaceType>
   static EIGEN_DEVICE_FUNC
-  void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, CoeffVectorType&, bool extractQ)
+  void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, CoeffVectorType&, WorkSpaceType&, bool extractQ)
   {
     diag(0,0) = numext::real(mat(0,0));
     if(extractQ)