From edd4dbe58857d77c239f47e3867ae2aa5ef4c8d4 Mon Sep 17 00:00:00 2001 From: Simon Mataigne <61558943+smataigne@users.noreply.github.com> Date: Thu, 18 Aug 2022 13:29:07 -0400 Subject: [PATCH] Revert "possible cleanups (#66)" (#67) This reverts commit 338103351bd235f2d1b76117f98920182551c96b. --- src/hessenberg.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/hessenberg.jl b/src/hessenberg.jl index 7c4409a..e586292 100644 --- a/src/hessenberg.jl +++ b/src/hessenberg.jl @@ -13,26 +13,26 @@ end LA.hessenberg(A::SkewHermitian)=hessenberg!(copyeigtype(A)) -# todo: add a comment explaining what this function returns + @views function householder!(x::AbstractVector{T},n::Integer) where {T} if n == 1 && T <:Real - return T(0), real(x[1]) # no final 1x1 reflection in the real case + return convert(T, 0), x[1] end - xnorm = norm(x[2:end]) + xnorm = (n > 1 ? norm(x[2:end]) : zero(real(x[1]))) alpha = x[1] - if !iszero(xnorm) || - (n == 1 && !iszero(alpha)) # sign reflection for 1x1 complex case - beta = (real(alpha) > 0 ? -1 : +1) * hypot(alpha,xnorm) + if !iszero(xnorm) || (n == 1 && !iszero(alpha)) + beta = (real(alpha) > 0 ? -1 : +1) * hypot(abs(alpha),xnorm) tau = 1 - alpha / beta alpha = 1 / (alpha - beta) x[1] = 1 x[2:n] .*= alpha + alpha = T(beta) else tau = T(0) - x .= 0 - beta = real(T)(0) + x .= zeros(T, n) + alpha = T(0) end - return tau, beta + return tau, alpha end @views function ger2!(tau::Number , v::StridedVector{T} , s::StridedVector{T}, @@ -101,7 +101,7 @@ end #Generate elementary reflector H(i) to annihilate A(i+2:n,i) stau,alpha = householder!(A[i+1:n,i],n-i) - E[i] = alpha + E[i] = real(alpha) mul!(W[i+1:n,i], A[i+1:n,i+1:n], A[i+1:n,i], 1, 0) if i>1 @@ -144,7 +144,7 @@ end nb = setnb(n) A = S.data - E = similar(A, real(T), n - 1) + E = similar(A, n - 1) tau = similar(A, n - 1) W = similar(A, n, nb) update = similar(A, n - nb, n - nb)