2x2 hermitian `eigen` method prone to under/overflow

[mikmoore]

julia> eigen(Hermitian(SA[2 1; 1 1] * 1f-30))
Eigen{Float32, Float32, SMatrix{2, 2, Float32, 4}, SVector{2, Float32}}
values:
2-element SVector{2, Float32} with indices SOneTo(2):
 3.8196598f-31
 2.618034f-30
vectors:
2×2 SMatrix{2, 2, Float32, 4} with indices SOneTo(2)×SOneTo(2):
 -Inf  Inf
  Inf  Inf

julia> eigen(Hermitian(SA[2 1; 1 1] * 1f+30))
Eigen{Float32, Float32, SMatrix{2, 2, Float32, 4}, SVector{2, Float32}}
values:
2-element SVector{2, Float32} with indices SOneTo(2):
 3.819661f29
 2.618034f30
vectors:
2×2 SMatrix{2, 2, Float32, 4} with indices SOneTo(2)×SOneTo(2):
 -0.0  0.0
  0.0  0.0

The relevant method appears to be function _eig(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}.

It looks like the issue is the lines n1 = sqrt(v11' * v11 + a21 * a21') and n2 = sqrt(v21' * v21 + a21 * a21'), which can under/overflow for small/large values. It appears that replacing them with something like the line tmp = norm(SVector(diag_avg_diff, a21)) (from further up the function) would resolve the issue. Should be very quick.

[martinholters]

Should that rather use hypot?

[mikmoore]

It definitely could (assuming it works for Complex), although hypot does extra work to be correctly rounded which is probably overkill here.

Another option is to compute the max(abs(arg1), abs(arg2)) and divide both arguments by that before proceeding with the calculation as written (although replacing x'*x with abs2(x) would also be an improvement). Without looking, this is probably pretty close to what norm does here.