From 9699cb0e61abd3e7978300ddeb7439b94c87c8ea Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 25 Jul 2023 14:20:03 +0530
Subject: [PATCH 1/2] Faster eigvecs for SymTridiagonal

---
 src/eigen.jl | 46 ++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 44 insertions(+), 2 deletions(-)

diff --git a/src/eigen.jl b/src/eigen.jl
index 14a36dd..4f76693 100644
--- a/src/eigen.jl
+++ b/src/eigen.jl
@@ -79,9 +79,9 @@ function _eigvecs_toeplitz(T)
     c1 = T[1,2]  # superdiagonal
     prefactors = _eigvec_prefactors(T, cm1, c1)
     for q in axes(M,2)
-        qrev = n+1-q # match the default eigenvalue sorting
         for j in 1:cld(n,2)
-            M[j, q] = prefactors[j] * sinpi(j*qrev/(n+1))
+            jphase = 2isodd(j) - 1
+            M[j, q] = prefactors[j] * jphase * sinpi(j * q/(n+1))
         end
         phase = iseven(n+q) ? 1 : -1
         for j in cld(n,2)+1:n
@@ -92,6 +92,48 @@ function _eigvecs_toeplitz(T)
     return M
 end
 
+function _eigvecs_toeplitz(T::Union{SymTridiagonal, Symmetric{<:Any,<:Tridiagonal}})
+    require_one_based_indexing(T)
+    n = checksquare(T)
+    M = Matrix{_eigvec_eltype(T)}(undef, n, n)
+    n == 0 && return M
+    n == 1 && return fill!(M, oneunit(eltype(M)))
+    for q in 1:cld(n,2)
+        for j in 1:q
+            jphase = 2isodd(j) - 1
+            M[j, q] = jphase * sinpi(j * q/(n+1))
+        end
+    end
+    for q in 1:cld(n,2)
+        for j in q+1:cld(n,2)
+            qphase = 2isodd(q) - 1
+            jphase = 2isodd(j) - 1
+            phase = qphase * jphase
+            M[j, q] = phase * M[q, j]
+        end
+    end
+    for q in 1:cld(n,2)
+        phase = iseven(n+q) ? 1 : -1
+        for j in cld(n,2)+1:n
+            M[j, q] = phase * M[n+1-j,q]
+        end
+    end
+    for q in cld(n,2)+1:n
+        for j in 1:cld(n,2)
+            qphase = 2isodd(q) - 1
+            jphase = 2isodd(j) - 1
+            phase = qphase * jphase
+            M[j, q] = phase * M[q, j]
+        end
+        phase = iseven(n+q) ? 1 : -1
+        for j in cld(n,2)+1:n
+            M[j, q] = phase * M[n+1-j,q]
+        end
+    end
+    _normalizecols!(M, T)
+    return M
+end
+
 function _eigvecs(A)
     n = size(A,1)
     if n <= 2 # repeated roots possible

From 977ad14b6783df9e81add00f806b630f41037122 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 25 Jul 2023 15:08:35 +0530
Subject: [PATCH 2/2] Test for more sizes

---
 src/eigen.jl     | 7 ++++---
 test/runtests.jl | 2 +-
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/eigen.jl b/src/eigen.jl
index 4f76693..739dba6 100644
--- a/src/eigen.jl
+++ b/src/eigen.jl
@@ -44,6 +44,10 @@ function _eigvals_toeplitz(T)
     return vals
 end
 
+_eigvec_eltype(A::Union{SymTridiagonal,
+                    Symmetric{<:Any,<:Tridiagonal}}) = float(eltype(A))
+_eigvec_eltype(A) = complex(float(eltype(A)))
+
 _eigvec_prefactor(A, cm1, c1, m) = sqrt(complex(cm1/c1))^m
 _eigvec_prefactor(A::Union{SymTridiagonal, Symmetric{<:Any, <:Tridiagonal}}, cm1, c1, m) = oneunit(_eigvec_eltype(A))
 
@@ -54,9 +58,6 @@ end
 _eigvec_prefactors(A::Union{SymTridiagonal, Symmetric{<:Any, <:Tridiagonal}}, cm1, c1) =
     Fill(_eigvec_prefactor(A, cm1, c1, 1), size(A,1))
 
-_eigvec_eltype(A::SymTridiagonal) = float(eltype(A))
-_eigvec_eltype(A) = complex(float(eltype(A)))
-
 @static if !isdefined(Base, :eachcol)
     eachcol(A) = (view(A,:,i) for i in axes(A,2))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 27ec70f..fd8e49d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -584,7 +584,7 @@ end
 @testset "eigen" begin
     sortby = x -> (real(x), imag(x))
     @testset "Tridiagonal Toeplitz" begin
-        _sizes = (1, 2, 6, 10)
+        _sizes = (1, 2, 5, 6, 10, 15)
         sizes = VERSION >= v"1.6" ? (0, _sizes...) : _sizes
         @testset for n in sizes
             @testset "Tridiagonal" begin