Fix size of Householder (#101)

Also improve error messages and extend Base's convert instead of defining a new one

Author: andreasnoack

src/householder.jl

@@ -10,8 +10,8 @@ struct HouseholderBlock{T,S<:StridedMatrix,U<:StridedMatrix}
     T::UpperTriangular{T,U}
 end
 
-size(H::Householder) = (length(H.v), length(H.v))
-size(H::Householder, i::Integer) = i <= 2 ? length(H.v) : 1
+size(H::Householder) = (length(H.v) + 1, length(H.v) + 1)
+size(H::Householder, i::Integer) = i <= 2 ? length(H.v) + 1 : 1
 
 eltype(H::Householder{T})      where T = T
 eltype(H::HouseholderBlock{T}) where T = T
@@ -21,7 +21,7 @@ adjoint(H::HouseholderBlock{T}) where {T} = Adjoint{T,typeof(H)}(H)
 
 function lmul!(H::Householder, A::StridedMatrix)
     m, n = size(A)
-    length(H.v) == m - 1 || throw(DimensionMismatch(""))
+    length(H.v) == m - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but first dimension of matrix is $(size(A, 1))"))
     v = view(H.v, 1:m - 1)
     τ = H.τ
     for j = 1:n
@@ -37,7 +37,7 @@ end
 
 function rmul!(A::StridedMatrix, H::Householder)
     m, n = size(A)
-    length(H.v) == n - 1 || throw(DimensionMismatch(""))
+    length(H.v) == n - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but second dimension of matrix is $(size(A, 2))"))
     v = view(H.v, :)
     τ = H.τ
     a1 = view(A, :, 1)
@@ -52,7 +52,7 @@ end
 function lmul!(adjH::Adjoint{<:Any,<:Householder}, A::StridedMatrix)
     H = parent(adjH)
     m, n = size(A)
-    length(H.v) == m - 1 || throw(DimensionMismatch(""))
+    length(H.v) == m - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but first dimension of matrix is $(size(A, 1))"))
     v = view(H.v, 1:m - 1)
     τ = H.τ
     for j = 1:n
@@ -142,5 +142,5 @@ end
 (*)(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}) where {T} =
         lmul!(adjH, copy(A), similar(A, (min(size(parent(adjH).V)...), size(A, 2))))
 
-convert(::Type{Matrix}, H::Householder{T}) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1)))
-convert(::Type{Matrix{T}}, H::Householder{T}) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1)))
+Base.convert(::Type{Matrix}, H::Householder{T}) where {T} = convert(Matrix{T}, H)
+Base.convert(::Type{Matrix{T}}, H::Householder) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1)))

test/eigengeneral.jl

@@ -144,4 +144,16 @@ end
     @test sort(imag(vals)) ≈ sort(imag(λs)) atol=1e-25
 end
 
+@testset "_hessenberg! and Hessenberg" begin
+    n = 10
+    A = randn(n, n)
+    HF = GenericLinearAlgebra._hessenberg!(copy(A))
+    for i in 1:length(HF.τ)
+        HM = convert(Matrix, HF.τ[i])
+        A[(i + 1):end, :] = HM * A[(i + 1):end, :]
+        A[:, (i + 1):end] = A[:, (i + 1):end] * HM'
+    end
+    @test tril(A, -2) ≈ zeros(n, n) atol = 1e-14
+end
+
 end