Finish eig()

* We currently support Hermitian eigenvalue decomposition, even for
  those not wrapped in `Hermition{}`
* Non-hermitian static matrices now error on eig() (there is a run-time
  test)
* Some matrix_multiply `@inbounds`

Author: Andy Ferris

src/eigen.jl

@@ -3,27 +3,26 @@
 end
 
 @inline function eig{T, SM <: StaticMatrix}(A::Base.LinAlg.HermOrSym{T,SM}; permute::Bool=true, scale::Bool=true)
-    _eig(Size(A), A, permute, scale)
+    _eig(Size(SM), A, permute, scale)
 end
 
 @inline function _eig(s::Size, A::StaticMatrix, permute, scale)
-    # TODO This is not type stable: both this fast branch, and the implementation of `Base.eigfact`.
-    # Decision needed in how to proceed. See e.g. JuliaLang/julia#12304
+    # Only cover the hermitian branch, for now ast least
+    # This also solves some type-stability issues such as arise in Base
     if ishermitian(A)
        return _eig(s, Hermitian(A), permute, scale)
+    else
+       error("Only hermitian matrices are diagonalizable by *StaticArrays*. Non-Hermitian matrices should be converted to `Array` first.")
     end
-    eigen = eigfact(Array(A); permute=permute, scale=scale)
-    return (Size(Size(typeof(A))[1])(eigen.values), s(eigen.vectors)) # Return a SizedArray
 end
 
-
 @inline function _eig{T<:Real}(s::Size, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale)
     eigen = eigfact(Hermitian(Array(parent(A))); permute=permute, scale=scale)
     return (s(eigen.values), s(eigen.vectors)) # Return a SizedArray
 end
 
 
-@inline function _eig(::Size{(1,1)}, A, permute, scale)
+@inline function _eig{T<:Real}(::Size{(1,1)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale)
     @inbounds return (SVector{1,T}((A[1],)), eye(SMatrix{1,1,T}))
 end
 
@@ -148,66 +147,3 @@ end
         end
     end
 end
-
-
-
-
-# TODO: the non-symmetric case: type stable version (real -> real) since it is more useful to us!
-
-#=
-@generated function eig{T<:Real, SM <: StaticMatrix}(A::Base.LinAlg.RealHermSymComplexHerm{T,SM}; permute::Bool=true, scale::Bool=true)
-    if size(SM) == (1,1)
-        return quote
-            $(Expr(:meta, :inline))
-            @inbounds return (SVector{1,T}((A[1],)), eye(SMatrix{1,1,T}))
-        end
-    elseif size(SM) == (2,2)
-        return quote
-            $(Expr(:meta, :inline))
-            a = A.data
-
-            if m2.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.
-                vals = SVector(t_half - tmp, t_half + tmp)
-f
-                @inbounds if A[3] == 0
-                    @inbounds if A[3] == 0
-                        vecs = eye(SMatrix{2,2,T})
-                    else
-                        @inbounds vecs = @SMatrix [ A[3]          A[3]         ;
-                                                    vals[1]-A[1]  vals[2]-A[1] ]
-                    end
-                else
-                    @inbounds v11 = vals[1]-A[4]
-                    @inbounds n1 = sqrt(v11'*v11 + A[2]'*A[2])
-                    v11 = v11 / n1
-                    @inbounds v12 = A[2] / n1
-
-                    @inbounds v21 = vals[2]-A[4]
-                    @inbounds n2 = sqrt(v21'*v21 + A[2]'*A[2])
-                    v21 = v21 / n2
-                    @inbounds v22 = A[2] / n2
-
-                    vecs = @SMatrix [ v11  v21 ;
-                                      v12  v22 ]
-                end
-                return (vals,vecs)
-            else
-
-            end
-        end
-    elseif size(SM) == (3,3)
-        error("not implemented")
-    else
-        return quote
-            $(Expr(:meta, :inline))
-            eigen = eigfact(A)
-            return (eigen.values, eigen.vectors)
-        end
-    end
-end
-=#

src/matrix_multiply.jl

@@ -489,7 +489,7 @@ end
         @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
     end
 end
-7
+
 
 # TODO aliasing problems if c === b?
 @generated function A_mul_B!{T1,T2,T3}(c::StaticVector{T1}, A::StaticMatrix{T2}, b::StaticVector{T3})
@@ -767,11 +767,11 @@ end
 
     exprs = [:(C[$(sub2ind(s, k1, k2))] = $(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sA[1], k2 = 1:sB[2]]
 
-    return Expr(:block,
-        Expr(:meta,:inline),
-        vect_exprs...,
-        exprs...
-    )
+    return quote
+        Expr(:meta,:inline)
+        @inbounds $(Expr(:block, vect_exprs...))
+        @inbounds $(Expr(:block, exprs...))
+    end
 end
 
 #function A_mul_B_blas(a, b, c, A, B)

test/eigen.jl

@@ -0,0 +1,42 @@
+@testset "Eigenvalue decomposition" begin
+    @testset "1×1" begin
+        m = @SMatrix [2.0]
+        (vals, vecs) = eig(m)
+        @test vals === SVector(2.0)
+        @test vecs === SMatrix{1,1}(1.0)
+
+        (vals, vecs) = eig(Symmetric(m))
+        @test vals === SVector(2.0)
+        @test vecs === SMatrix{1,1}(1.0)
+    end
+
+    @testset "2×2" for i = 1:100
+        m_a = randn(2,2)
+        m_a = m_a*m_a'
+        m = SMatrix{2,2}(m_a)
+
+        (vals_a, vecs_a) = eig(m)
+        (vals, vecs) = eig(m)
+        @test vals::SVector ≈ vals_a
+        @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+
+        (vals, vecs) = eig(Symmetric(m))
+        @test vals::SVector ≈ vals_a
+        @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+    end
+
+    @testset "3×3" for i = 1:100
+        m_a = randn(3,3)
+        m_a = m_a*m_a'
+        m = SMatrix{3,3}(m_a)
+
+        (vals_a, vecs_a) = eig(m)
+        (vals, vecs) = eig(m)
+        @test vals::SVector ≈ vals_a
+        @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+
+        (vals, vecs) = eig(Symmetric(m))
+        @test vals::SVector ≈ vals_a
+        @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+    end
+end

test/matrix_multiply.jl

@@ -1,4 +1,4 @@
-@testset "Matrix multiplication" begin
+@testset "Matrix multiplication" begin 
     @testset "Matrix-vector" begin
         m = @SMatrix [1 2; 3 4]
         v = @SVector [1, 2]
@@ -163,46 +163,46 @@
         @test a2::MArray{(2,2),Int,2,4} == @MArray [10 13; 22 29]
 
         # Alternative builtin method used for n > 8
-        m_array = rand(1:10, 10, 10)
-        n_array = rand(1:10, 10, 10)
-        a_array = m_array*n_array
+        m_array_2 = rand(1:10, 10, 10)
+        n_array_2 = rand(1:10, 10, 10)
+        a_array_2 = m_array_2*n_array_2
 
-        m = MMatrix{10,10}(m_array)
-        n = MMatrix{10,10}(n_array)
-        a = MMatrix{10,10,Int}()
-        A_mul_B!(a, m, n)
-        @test a == a_array
+        m_2 = MMatrix{10,10}(m_array_2)
+        n_2 = MMatrix{10,10}(n_array_2)
+        a_2 = MMatrix{10,10,Int}()
+        A_mul_B!(a_2, m_2, n_2)
+        @test a_2 == a_array_2
 
         # BLAS used for n > 14
-        m_array = randn(4, 4)
-        n_array = randn(4, 4)
-        a_array = m_array*n_array
-
-        m = MMatrix{4,4}(m_array)
-        n = MMatrix{4,4}(n_array)
-        a = MMatrix{4,4,Float64}()
-        A_mul_B!(a, m, n)
-        @test a ≈ a_array
-
-        m_array = randn(10, 10)
-        n_array = randn(10, 10)
-        a_array = m_array*n_array
-
-        m = MMatrix{10,10}(m_array)
-        n = MMatrix{10,10}(n_array)
-        a = MMatrix{10,10,Float64}()
-        A_mul_B!(a, m, n)
-        @test a ≈ a_array
-
-        m_array = rand(1:10, 16, 16)
-        n_array = rand(1:10, 16, 16)
-        a_array = m_array*n_array
-
-        m = MMatrix{16,16}(m_array)
-        n = MMatrix{16,16}(n_array)
-        a = MMatrix{16,16,Int}()
-        A_mul_B!(a, m, n)
-        @test a ≈ a_array
+        m_array_3 = randn(4, 4)
+        n_array_3 = randn(4, 4)
+        a_array_3 = m_array_3*n_array_3
+
+        m_3 = MMatrix{4,4}(m_array_3)
+        n_3 = MMatrix{4,4}(n_array_3)
+        a_3 = MMatrix{4,4,Float64}()
+        A_mul_B!(a_3, m_3, n_3)
+        @test a_3 ≈ a_array_3
+
+        m_array_4 = randn(10, 10)
+        n_array_4 = randn(10, 10)
+        a_array_4 = m_array_4*n_array_4
+
+        m_4 = MMatrix{10,10}(m_array_4)
+        n_4 = MMatrix{10,10}(n_array_4)
+        a_4 = MMatrix{10,10,Float64}()
+        A_mul_B!(a_4, m_4, n_4)
+        @test a_4 ≈ a_array_4
+
+        m_array_5 = rand(1:10, 16, 16)
+        n_array_5 = rand(1:10, 16, 16)
+        a_array_5 = m_array_5*n_array_5
+
+        m_5 = MMatrix{16,16}(m_array_5)
+        n_5 = MMatrix{16,16}(n_array_5)
+        a_5 = MMatrix{16,16,Int}()
+        A_mul_B!(a_5, m_5, n_5)
+        @test a_5 ≈ a_array_5
 
         # Float64
         vf = @SVector [2.0, 4.0]

test/runtests.jl

@@ -19,5 +19,6 @@ using Base.Test
     include("linalg.jl")
     include("matrix_multiply.jl")
     include("solve.jl")
+    include("eigen.jl")
     include("deque.jl")
 end