Merge pull request #180 from kshyatt/ksh/eigenagain

Let eigen work for "big" Herm/symm matrices

Author: c42f

src/eigen.jl

@@ -12,6 +12,7 @@ end
 end
 
 @inline _eigvals(::Size{(1,1)}, a, permute, scale) = @inbounds return SVector(real(a.data[1]))
+@inline _eigvals(::Size{(1, 1)}, a::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} = @inbounds return SVector(real(parent(a).data[1]))
 
 @inline function _eigvals(::Size{(2,2)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
     a = A.data
@@ -102,6 +103,10 @@ end
     return SVector(eig1, eig2, eig3)
 end
 
+@inline function _eigvals(s::Size, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
+    vals = eigvals(Hermitian(Array(parent(A))))
+    return SVector{s[1], T}(vals)
+end
 
 
 
@@ -124,8 +129,8 @@ end
 end
 
 @inline function _eig(s::Size, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
-    eigen = eigfact(Hermitian(Array(parent(A))); permute=permute, scale=scale)
-    return (s(eigen.values), s(eigen.vectors)) # Return a SizedArray
+    eigen = eigfact(Hermitian(Array(parent(A))))
+    return (SVector{s[1], T}(eigen.values), SMatrix{s[1], s[2], T}(eigen.vectors)) # Return a SizedArray
 end
 
 

test/eigen.jl

@@ -141,4 +141,33 @@
         @test vecs*diagm(vals)*vecs' ≈ m
         @test eigvals(m) ≈ vals
     end
+    
+    @testset "4×4" for i = 1:100
+        m_a = randn(4,4)
+        m_a = m_a*m_a'
+        m = SMatrix{4,4}(m_a)
+
+        (vals_a, vecs_a) = eig(m_a)
+        (vals, vecs) = eig(m)
+        @test vals::SVector ≈ vals_a
+        @test eigvals(m) ≈ vals
+        @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+
+        (vals, vecs) = eig(Symmetric(m))
+        @test vals::SVector ≈ vals_a
+        @test eigvals(m) ≈ vals
+        @test eigvals(Hermitian(m)) ≈ vals
+        @test eigvals(Hermitian(m, :L)) ≈ vals
+        @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+       
+        (vals, vecs) = eig(Symmetric(m, :L))
+        @test vals::SVector ≈ vals_a
+        m_d = randn(SVector{4}); m = diagm(m_d)
+        (vals, vecs) = eig(Hermitian(m))
+        @test vals::SVector ≈ sort(m_d)
+        (vals, vecs) = eig(Hermitian(m, :L))
+        @test vals::SVector ≈ sort(m_d)
+        @test eigvals(m) ≈ sort(m_d)
+        @test eigvals(Hermitian(m)) ≈ sort(m_d)
+    end
 end