Added eigvals, tests

Andy Ferris · Andy Ferris · commit 261939b1f39c · 2016-12-02T10:38:57.000+10:00
diff --git a/src/StaticArrays.jl b/src/StaticArrays.jl
@@ -6,8 +6,9 @@ import Base: @pure, @propagate_inbounds, getindex, setindex!, size, similar,
              length, convert, promote_op, map, map!, reduce, mapreduce,
              broadcast, broadcast!, conj, transpose, ctranspose, hcat, vcat,
              ones, zeros, eye, one, cross, vecdot, reshape, fill, fill!, det,
-             inv, eig, trace, vecnorm, norm, dot, diagm, sum, prod, count, any,
-             all, sumabs, sumabs2, minimum, maximum, extrema, mean, copy
+             inv, eig, eigvals, trace, vecnorm, norm, dot, diagm, sum, prod,
+             count, any, all, sumabs, sumabs2, minimum, maximum, extrema, mean,
+             copy
 
 export StaticScalar, StaticArray, StaticVector, StaticMatrix
 export Scalar, SArray, SVector, SMatrix
diff --git a/src/eigen.jl b/src/eigen.jl
@@ -1,3 +1,107 @@
+
+@inline eigvals{T<:Real,SA<:StaticArray}(a::Base.LinAlg.RealHermSymComplexHerm{T,SA},; permute::Bool=true, scale::Bool=true) = _eigvals(Size(SA), a, permute, scale)
+@inline function eigvals(a::StaticArray; permute::Bool=true, scale::Bool=true)
+    if ishermitian(a)
+        _eigvals(Size(a), Hermitian(a), permute, scale)
+    else
+        error("Only hermitian matrices are diagonalizable by *StaticArrays*. Non-Hermitian matrices should be converted to `Array` first.")
+    end
+end
+
+@inline _eigvals(::Size{(1,1)}, a, permute, scale) = @inbounds return SVector(real(a.data[1]))
+
+@inline function _eigvals{T<:Real}(::Size{(2,2)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale)
+    a = A.data
+
+    if A.uplo == 'U'
+        @inbounds t_half = real(a[1] + a[4])/2
+        @inbounds d = real(a[1]*a[4] - a[3]'*a[3]) # Should be real
+
+        tmp2 = t_half*t_half - d
+        tmp2 < 0 ? tmp = zero(tmp2) : tmp = sqrt(tmp2) # Numerically stable for identity matrices, etc.
+        return SVector(t_half - tmp, t_half + tmp)
+    else
+        @inbounds t_half = real(a[1] + a[4])/2
+        @inbounds d = real(a[1]*a[4] - a[2]'*a[2]) # Should be real
+
+        tmp2 = t_half*t_half - d
+        tmp2 < 0 ? tmp = zero(tmp2) : tmp = sqrt(tmp2) # Numerically stable for identity matrices, etc.
+        return SVector(t_half - tmp, t_half + tmp)
+    end
+end
+
+@inline function _eigvals{T<:Real}(::Size{(3,3)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale)
+    S = typeof((one(T)*zero(T) + zero(T))/one(T))
+    Sreal = real(S)
+
+    @inbounds a11 = convert(Sreal, A.data[1])
+    @inbounds a22 = convert(Sreal, A.data[5])
+    @inbounds a33 = convert(Sreal, A.data[9])
+    if A.uplo == 'U'
+        @inbounds a12 = convert(S, A.data[4])
+        @inbounds a13 = convert(S, A.data[7])
+        @inbounds a23 = convert(S, A.data[8])
+    else
+        @inbounds a12 = conj(convert(S, A.data[2]))
+        @inbounds a13 = conj(convert(S, A.data[3]))
+        @inbounds a23 = conj(convert(S, A.data[6]))
+    end
+
+    p1 = abs2(a12) + abs2(a13) + abs2(a23)
+    if (p1 == 0)
+        # Matrix is diagonal
+        if a11 < a22
+            if a22 < a33
+                return SVector(a11, a22, a33)
+            elseif a33 < a11
+                return SVector(a33, a11, a22)
+            else
+                return SVector(a11, a33, a22)
+            end
+        else #a22 < a11
+            if a11 < a33
+                return SVector(a22, a11, a33)
+            elseif a33 < a22
+                return SVector(a33, a22, a11)
+            else
+                return SVector(a22, a33, a11)
+            end
+        end
+    end
+
+    q = (a11 + a22 + a33) / 3
+    p2 = abs2(a11 - q) + abs2(a22 - q) + abs2(a33 - q) + 2 * p1
+    p = sqrt(p2 / 6)
+    invp = inv(p)
+    b11 = (a11 - q) * invp
+    b22 = (a22 - q) * invp
+    b33 = (a33 - q) * invp
+    b12 = a12 * invp
+    b13 = a13 * invp
+    b23 = a23 * invp
+    B = SMatrix{3,3,S}((b11, conj(b12), conj(b13), b12, b22, conj(b23), b13, b23, b33))
+    r = real(det(B)) / 2
+
+    # In exact arithmetic for a symmetric matrix  -1 <= r <= 1
+    # but computation error can leave it slightly outside this range.
+    if (r <= -1)
+        phi = Sreal(pi) / 3
+    elseif (r >= 1)
+        phi = zero(Sreal)
+    else
+        phi = acos(r) / 3
+    end
+
+    eig3 = q + 2 * p * cos(phi)
+    eig1 = q + 2 * p * cos(phi + (2*Sreal(pi)/3))
+    eig2 = 3 * q - eig1 - eig3     # since trace(A) = eig1 + eig2 + eig3
+
+    return SVector(eig1, eig2, eig3)
+end
+
+
+
+
 @inline function eig(A::StaticMatrix; permute::Bool=true, scale::Bool=true)
     _eig(Size(A), A, permute, scale)
 end
@@ -106,8 +210,27 @@ end
     p1 = abs2(a12) + abs2(a13) + abs2(a23)
     if (p1 == 0)
         # Matrix is diagonal
-        # TODO need to sort the eigenvalues
-        return (SVector(a11, a22, a33), eye(SMatrix{3,3,S}))
+        v1 = SVector(one(S),  zero(S), zero(S))
+        v2 = SVector(zero(S), one(S),  zero(S))
+        v3 = SVector(zero(S), zero(S), one(S) )
+
+        if a11 < a22
+            if a22 < a33
+                return (SVector(a11, a22, a33), hcat(v1,v2,v3))
+            elseif a33 < a11
+                return (SVector(a33, a11, a22), hcat(v3,v1,v2))
+            else
+                return (SVector(a11, a33, a22), hcat(v1,v3,v2))
+            end
+        else #a22 < a11
+            if a11 < a33
+                return (SVector(a22, a11, a33), hcat(v2,v1,v3))
+            elseif a33 < a22
+                return (SVector(a33, a22, a11), hcat(v3,v2,v1))
+            else
+                return (SVector(a22, a33, a11), hcat(v2,v3,v1))
+            end
+        end
     end
 
     q = (a11 + a22 + a33) / 3
@@ -251,6 +374,8 @@ end
 #=
 Boost Software License - Version 1.0 - August 17th, 2003
 
+
+
 Permission is hereby granted, free of charge, to any person or organization
 obtaining a copy of the software and accompanying documentation covered by
 this license (the "Software") to use, reproduce, display, distribute,
diff --git a/test/eigen.jl b/test/eigen.jl
@@ -3,6 +3,7 @@
         m = @SMatrix [2.0]
         (vals, vecs) = eig(m)
         @test vals === SVector(2.0)
+        @test eigvals(m) === vals
         @test vecs === SMatrix{1,1}(1.0)
 
         (vals, vecs) = eig(Symmetric(m))
@@ -18,10 +19,12 @@
         (vals_a, vecs_a) = eig(m)
         (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 (vecs*diagm(vals)*vecs')::SMatrix ≈ m
     end
 
@@ -33,72 +36,81 @@
         (vals_a, vecs_a) = eig(m)
         (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 (vecs*diagm(vals)*vecs')::SMatrix ≈ m
     end
 
     @testset "3x3 degenerate cases" begin
         # Rank 1
         v = randn(SVector{3,Float64})
         m = v*v'
-        eigvals, eigvecs = eig(m)::Tuple{SVector,SMatrix}
+        vals, vecs = eig(m)::Tuple{SVector,SMatrix}
 
-        @test isapprox(eigvecs'*eigvecs, eye(SMatrix{3,3,Float64}); atol = 1e-4) # This algorithm isn't super accurate
-        @test eigvals ≈ SVector(0.0, 0.0, sumabs2(v))
+        @test vecs'*vecs ≈ eye(SMatrix{3,3,Float64})
+        @test vals ≈ SVector(0.0, 0.0, sumabs2(v))
+        @test eigvals(m) ≈ vals
 
         # Rank 2
         v2 = randn(SVector{3,Float64})
         v2 -= dot(v,v2)*v/sumabs2(v)
         m += v2*v2'
-        eigvals, eigvecs = eig(m)::Tuple{SVector,SMatrix}
+        vals, vecs = eig(m)::Tuple{SVector,SMatrix}
 
-        @test isapprox(eigvecs'*eigvecs, eye(SMatrix{3,3,Float64}); atol = 1e-4)
+        @test vecs'*vecs ≈ eye(SMatrix{3,3,Float64})
         if sumabs2(v) < sumabs2(v2)
-            @test eigvals ≈ SVector(0.0, sumabs2(v), sumabs2(v2))
+            @test vals ≈ SVector(0.0, sumabs2(v), sumabs2(v2))
         else
-            @test eigvals ≈ SVector(0.0, sumabs2(v2), sumabs2(v))
+            @test vals ≈ SVector(0.0, sumabs2(v2), sumabs2(v))
         end
+        @test eigvals(m) ≈ vals
 
         # Degeneracy (2 large)
         m = -99*(v*v')/sumabs2(v) + 100*eye(SMatrix{3,3,Float64})
-        eigvals, eigvecs = eig(m)::Tuple{SVector,SMatrix}
+        vals, vecs = eig(m)::Tuple{SVector,SMatrix}
 
-        @test isapprox(eigvecs'*eigvecs, eye(SMatrix{3,3,Float64}); atol = 1e-4)
-        @test eigvals ≈ SVector(1.0, 100.0, 100.0)
+        @test vecs'*vecs ≈ eye(SMatrix{3,3,Float64})
+        @test vals ≈ SVector(1.0, 100.0, 100.0)
+        @test eigvals(m) ≈ vals
 
         # Degeneracy (2 small)
         m = (v*v')/sumabs2(v) + 1e-2*eye(SMatrix{3,3,Float64})
-        eigvals, eigvecs = eig(m)::Tuple{SVector,SMatrix}
+        vals, vecs = eig(m)::Tuple{SVector,SMatrix}
 
-        @test isapprox(eigvecs'*eigvecs, eye(SMatrix{3,3,Float64}); atol = 1e-4)
-        @test eigvals ≈ SVector(1e-2, 1e-2, 1.01)
+        @test vecs'*vecs ≈ eye(SMatrix{3,3,Float64})
+        @test vals ≈ SVector(1e-2, 1e-2, 1.01)
+        @test eigvals(m) ≈ vals
 
         # Block diagonal
         m = @SMatrix [1.0 0.0 0.0;
                       0.0 1.0 1.0;
                       0.0 1.0 1.0]
-        eigvals, eigvecs = eig(m)::Tuple{SVector,SMatrix}
+        vals, vecs = eig(m)::Tuple{SVector,SMatrix}
 
-        @test eigvals ≈ [0.0, 1.0, 2.0]
-        @test eigvecs*diagm(eigvals)*eigvecs' ≈ m
+        @test vals ≈ [0.0, 1.0, 2.0]
+        @test vecs*diagm(vals)*vecs' ≈ m
+        @test eigvals(m) ≈ vals
 
         m = @SMatrix [1.0 0.0 1.0;
                       0.0 1.0 0.0;
                       1.0 0.0 1.0]
-        eigvals, eigvecs = eig(m)::Tuple{SVector,SMatrix}
+        vals, vecs = eig(m)::Tuple{SVector,SMatrix}
 
-        @test eigvals ≈ [0.0, 1.0, 2.0]
-        @test eigvecs*diagm(eigvals)*eigvecs' ≈ m
+        @test vals ≈ [0.0, 1.0, 2.0]
+        @test vecs*diagm(vals)*vecs' ≈ m
+        @test eigvals(m) ≈ vals
 
         m = @SMatrix [1.0 1.0 0.0;
                       1.0 1.0 0.0;
                       0.0 0.0 1.0]
-        eigvals, eigvecs = eig(m)::Tuple{SVector,SMatrix}
+        vals, vecs = eig(m)::Tuple{SVector,SMatrix}
 
-        @test eigvals ≈ [0.0, 1.0, 2.0]
-        @test eigvecs*diagm(eigvals)*eigvecs' ≈ m
+        @test vals ≈ [0.0, 1.0, 2.0]
+        @test vecs*diagm(vals)*vecs' ≈ m
+        @test eigvals(m) ≈ vals
     end
 end
