From 56006774fb88562417ca4ced2552a4ea4706b0c6 Mon Sep 17 00:00:00 2001
From: Aravindh Krishnamoorthy <aravindh.krishnamoorthy@fau.de>
Date: Fri, 7 Jul 2023 19:44:53 +0200
Subject: [PATCH] Memory optimised versions

---
 stdlib/LinearAlgebra/src/bunchkaufman.jl   |  2 +-
 stdlib/LinearAlgebra/src/symmetriceigen.jl | 46 ++++++++++++++++++----
 2 files changed, 39 insertions(+), 9 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/bunchkaufman.jl b/stdlib/LinearAlgebra/src/bunchkaufman.jl
index d1019a1a4ea5a..9f51836f92ca9 100644
--- a/stdlib/LinearAlgebra/src/bunchkaufman.jl
+++ b/stdlib/LinearAlgebra/src/bunchkaufman.jl
@@ -88,7 +88,7 @@ Base.iterate(S::BunchKaufman) = (S.D, Val(:UL))
 Base.iterate(S::BunchKaufman, ::Val{:UL}) = (S.uplo == 'L' ? S.L : S.U, Val(:p))
 Base.iterate(S::BunchKaufman, ::Val{:p}) = (S.p, Val(:done))
 Base.iterate(S::BunchKaufman, ::Val{:done}) = nothing
-
+copy(S::BunchKaufman) = BunchKaufman(copy(S.LD), copy(S.ipiv), S.uplo, S.symmetric, S.rook, S.info)
 
 """
     bunchkaufman!(A, rook::Bool=false; check = true) -> BunchKaufman
diff --git a/stdlib/LinearAlgebra/src/symmetriceigen.jl b/stdlib/LinearAlgebra/src/symmetriceigen.jl
index f92d0311e9c92..2982284707a6b 100644
--- a/stdlib/LinearAlgebra/src/symmetriceigen.jl
+++ b/stdlib/LinearAlgebra/src/symmetriceigen.jl
@@ -186,16 +186,31 @@ end
 
 # Bunch-Kaufmann (LDLT) based solution for generalized eigenvalues and eigenvectors
 function eigen(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
-    eigen!(copy(A), B; sortby)
+    eigen!(copy(A), copy(B); sortby)
 end
 function eigen!(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
-    D = B.D
-    M = (B.uplo == 'U') ? B.U : B.L ;
+   # Call sysconvf/_rook directly as getproperty creates a copy of B.LD.
+   if B.rook
+        LUD, od = LAPACK.syconvf_rook!(B.uplo, 'C', B.LD, B.ipiv)
+    else
+        LUD, od = LAPACK.syconv!(B.uplo, B.LD, B.ipiv)
+    end
+    if B.uplo == 'U'
+        M = UnitUpperTriangular(LUD)
+        du = od[2:end]
+        # Aliasing pointers of dl and du is not allowed for gtsv!
+        dl = B.symmetric ? copy(du) : conj.(copy(du))
+    else
+        M = UnitLowerTriangular(LUD)
+        dl = od[1:end-1]
+        # Aliasing pointers of dl and du is not allowed for gtsv!
+        du = B.symmetric ? copy(dl) : conj.(copy(dl))
+    end
     LAPACK.lapmt!(A, B.p, true)
     LAPACK.lapmr!(A, B.p, true)
     LAPACK.trtrs!(B.uplo, 'N', 'U', M.data, A)
     BLAS.trsm!('R', B.uplo , 'C', 'U', one(T), M.data, A)
-    LAPACK.gtsv!(D.dl, D.d, copy(D.du), A)
+    LAPACK.gtsv!(dl, diag(LUD), du, A)
     vals, vecs = eigen!(A; sortby)
     LAPACK.trtrs!(B.uplo, 'C', 'U', convert(typeof(vecs), M.data), vecs)
     LAPACK.lapmr!(vecs, invperm(B.p), true)
@@ -239,15 +254,30 @@ end
 
 # Bunch-Kaufmann (LDLT) based solution for generalized eigenvalues
 function eigvals(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
-    eigvals!(copy(A), B; sortby)
+    eigvals!(copy(A), copy(B); sortby)
 end
 function eigvals!(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
-    D = B.D
-    M = (B.uplo == 'U') ? B.U : B.L ;
+    # Call sysconvf/_rook directly as getproperty creates a copy of B.LD.
+    if B.rook
+        LUD, od = LAPACK.syconvf_rook!(B.uplo, 'C', B.LD, B.ipiv)
+    else
+        LUD, od = LAPACK.syconv!(B.uplo, B.LD, B.ipiv)
+    end
+    if B.uplo == 'U'
+        M = UnitUpperTriangular(LUD)
+        du = od[2:end]
+        # Aliasing pointers of dl and du is not allowed for gtsv!
+        dl = B.symmetric ? copy(du) : conj.(copy(du))
+    else
+        M = UnitLowerTriangular(LUD)
+        dl = od[1:end-1]
+        # Aliasing pointers of dl and du is not allowed for gtsv!
+        du = B.symmetric ? copy(dl) : conj.(copy(dl))
+    end
     LAPACK.lapmt!(A, B.p, true)
     LAPACK.lapmr!(A, B.p, true)
     LAPACK.trtrs!(B.uplo, 'N', 'U', M.data, A)
     BLAS.trsm!('R', B.uplo , 'C', 'U', one(T), M.data, A)
-    LAPACK.gtsv!(D.dl, D.d, copy(D.du), A)
+    LAPACK.gtsv!(dl, diag(LUD), du, A)
     return eigvals!(A; sortby)
 end