Handle non-hermitian matrices in eigen

Author: CarpeNecopinum

src/eigen.jl

@@ -124,13 +124,14 @@ end
 end
 
 
-@inline function _eig(s::Size, A::StaticMatrix, permute, scale)
-    # Only cover the hermitian branch, for now at least
-    # This also solves some type-stability issues that arise in Base
+@inline function _eig(s::Size, A::T, permute, scale) where {T <: StaticMatrix}
+    # For the non-hermitian branch, fall back to LinearAlgebra
     if ishermitian(A)
-       return _eig(s, Hermitian(A), permute, scale)
+        eivals, eivecs = _eig(s, Hermitian(A), permute, scale)
+        return eivals, eivecs
     else
-       error("Only hermitian matrices are diagonalizable by *StaticArrays*. Non-Hermitian matrices should be converted to `Array` first.")
+        eivals, eivecs = eigen(Array(A); permute = permute, scale = scale)
+        return SVector{s[1]}(eivals), SMatrix{s[1],s[2]}(eivecs)
     end
 end
 

test/eigen.jl

@@ -232,9 +232,6 @@ using StaticArrays, Test, LinearAlgebra
         @test vals::SVector ≈ sort(m_d)
         @test eigvals(m) ≈ sort(m_d)
         @test eigvals(Hermitian(m)) ≈ sort(m_d)
-
-        # not Hermitian
-        @test_throws Exception eigen(@SMatrix randn(4,4))
     end
 
     @testset "complex" begin
@@ -248,4 +245,45 @@ using StaticArrays, Test, LinearAlgebra
             @test V'*A*V ≈ diagm(Val(0) => D)
         end
     end
+
+    @testset "hermitian type stability" begin
+        for n=1:4
+            m = @SMatrix randn(n,n)
+            m += m'
+
+            @inferred eigen(Hermitian(m))
+            @inferred eigen(Symmetric(m))
+
+            mc = @SMatrix randn(ComplexF64, n, n)
+            @inferred eigen(Hermitian(mc + mc'))
+        end
+    end
+
+    @testset "non-hermitian 2d" begin
+        for n=1:5
+            angle = 2π * rand()
+            rot = @SMatrix [cos(angle) -sin(angle); sin(angle) cos(angle)]
+
+            vals, vecs = eigen(rot)
+
+            @test norm(vals[1]) ≈ 1.0
+            @test norm(vals[2]) ≈ 1.0
+
+            @test vecs[:,1] ≈ conj.(vecs[:,2])
+        end
+    end
+
+    @testset "non-hermitian 3d" begin
+        for n=1:5
+            angle = 2π * rand()
+            rot = @SMatrix [cos(angle) 0.0 -sin(angle); 0.0 1.0 0.0; sin(angle) 0.0 cos(angle)]
+
+            vals, vecs = eigen(rot)
+
+            @test norm(vals[1]) ≈ 1.0
+            @test norm(vals[2]) ≈ 1.0
+
+            @test vecs[:,1] ≈ conj.(vecs[:,2])
+        end
+    end
 end