(RFC) fix eigen in the diagonal 2x2-case

[dkarrasch]

Fixes #523. I also reduced the amount of @inbounds by tearing it to the outside. I think it should apply everywhere in the following block.

Some questions:

• Is T equal to eltype(A)?
• Is there a noticeable difference between TA(1), TA(0) and one(TA), zero(TA)?

TODO: add tests for the diagonal 2x2 case.

[andyferris]

Is there a noticeable difference between TA(1), TA(0) and one(TA), zero(TA)?

Yes, TA(1) is a constructor and might do whatever that type wants to do with a constructor, it is probably better to use convert(TA, 1) (see JuliaLang/julia#29350).

I'm guessing for Number this will tend to agree with one(TA). However, in general one returns a multiplicative identity, and some types like square, static matrices might define one as the identity matrix, but error out on convert(SMatrix{2,2,Float64,4), 1).

one and zero are possibly more likely to return constants than convert, though I suspect in most cases the compiler and LLVM would produce exactly the same code.

This code is for Real, so preferably use zero, one, and convert is fine - but at this stage I would recommend not to use the constructor.

[dkarrasch]

Thanks, @andyferris . I just checked Base, one found that
https://github.com/JuliaLang/julia/blob/70be5f25981918af138cce21bace331aa2d16945/base/number.jl#L297
so these don't make a difference, in contrast to the constructors. I'll replace them with one/zero. The other question is whether the T in A::LinearAlgebra.RealHermSymComplexHerm{T} always coincides with eltype(A)? If that's the case, then I could reduce the amount of types and variables floating around a bit.

[fredrikekre]

Unrelated to the bugfix, but the codepath for 'L' and 'U' is only different in 2 places, and I think this method could just be

@inline function _eig(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
    a = A.data
    TA = eltype(A)
    lower = A.uplo === 'L'

    @inbounds begin
        A11 = a[1]
        A22 = a[4]
        A2112 = lower ? a[2] : a[3]

        if A2112 == 0 # A is diagonal
            if A11 < A22
                vals = SVector(A11, A22)
                vecs = @SMatrix [one(TA) zero(TA);
                                 zero(TA) one(TA)]
            else # A22 <= A11
                vals = SVector(A22, A11)
                vecs = @SMatrix [zero(TA) one(TA);
                                 one(TA) zero(TA)]
            end
        else # A is not diagonal
            t_half = real(A11 + A22) / 2
            d = real(A11 * A22 - A2112' * A2112) # Should be real

            tmp2 = t_half * t_half - d
            tmp = tmp2 < 0 ? zero(tmp2) : sqrt(tmp2) # Numerically stable for identity matrices, etc.
            vals = SVector(t_half - tmp, t_half + tmp)

            v11 = vals[1] - A22
            n1 = sqrt(v11' * v11 + A2112' * A2112)
            v11 = v11 / n1
            v12 = (lower ? A2112 : A2112') / n1

            v21 = vals[2] - A22
            n2 = sqrt(v21' * v21 + A2112' * A2112)
            v21 = v21 / n2
            v22 = (lower ? A2112 : A2112') / n2

            vecs = @SMatrix [ v11  v21 ;
                              v12  v22 ]
        end
    end
    return (vals, vecs)
end

which reduced the code duplication quite a lot.

[dkarrasch]

Yes, I was aware of that. I wasn't sure about the original intention to have that code duplication, so I felt it must have been something beyond my understanding of compiler functionality (which is very low). I'll use @fredrikekre's version then.

[fredrikekre]

It was just a suggestion, since code duplication can often lead to bugs. I should perhaps have benchmarked first: StaticArrays master results in 7.8 ns and my suggestion 8.0 ns so perhaps worth keeping the big if-else-end block.

[dkarrasch]

Hm, on my machine (tested on a Symmetric(SMatrix)) I get 7.290 ns on master vs 7.532 ns on this PR vs. 7.624 ns with @fredrikekre's suggestion. So it seems that we do pay a price here, anyways. I can reduce runtime slightly to 7.512 ns by replacing a[*] == 0 by a[*] == zero(TA)!

[mschauer]

Is this ever needed? You can compare with 0, won't that work fine? For example no need to create a zero(DualNumber) if the type is a dual number.

[dkarrasch]

I guess it depends on what is cheaper: converting 0 to a comparable type or constructing the zero element. I just tested one versus the other, and for Float it doesn't make any difference (in contrast to my minor improvement above), but for Dual the comparison == zero(Dual) is indeed more expensive than == 0. I guess I'll revert the last commit then?

[fredrikekre]

I added some tests, we can refactor the code some other time, but it is good to get in the bugfix ASAP