Fix real matrix sqrt and log for edge cases

[sethaxen]

Fixes JuliaLang/LinearAlgebra.jl#825, also eliminates an unnecessary scalar sqrt tries to reduce potential numerical issues in _sqrt_real_2x2!.

The added tests fail on master. Test 1 fails unless we make this change to _sqrt_real_2x2!. Test 2 also fails unless we make this change to _sqrt_quasitriu_offdiag_block_2x2!.

Edit: Test 3 tests the case where NaNs would cause LAPACK.trsyl! to error.

[oscardssmith]

How does this affect the performance?

[sethaxen]

How does this affect the performance?

This version of _sqrt_real_2x2! takes about 5ns longer than on master (~35ns on master, ~40ns here). Using the same benchmark code and matrices as in #39973 (comment), sqrt_quasitriu takes about 16.5 μs on master and 17.2μs here.

[sethaxen]

With these changes, _sqrt_real_2x2! is ~33% slower. That being said, on a 20x20 matrix, sqrt_quasitriu is only ~4% slower. sqrt overflows and underflows much less with these changes, so I think this decrease in performance is worth it.

Updated Benchmarks

Setup:

using BenchmarkTools, LinearAlgebra
using Random; Random.seed!(1)
F = schur(randn(20, 20)^2)
x = F.T[1:2,1:2] # 2x2 block
y = zero(x)

Before:

julia> @btime $(LinearAlgebra._sqrt_real_2x2!)($y, $x);
  35.130 ns (0 allocations: 0 bytes)
 
julia> @btime $(LinearAlgebra.sqrt_quasitriu)($(F.T));
  16.880 μs (1 allocation: 3.25 KiB)

This PR:

julia> @btime $(LinearAlgebra._sqrt_real_2x2!)($y, $x);
  47.072 ns (0 allocations: 0 bytes)

julia> @btime $(LinearAlgebra.sqrt_quasitriu)($(F.T));
  17.599 μs (1 allocation: 3.25 KiB)

[sethaxen]

The 2x2 diagonal block computation of log suffered from the same overflow/underflow, so I also adopted a similar approach for that.

[sethaxen]

With the recent changes, _sqrt_real_2x2! is almost as fast as the original:

julia> @btime $(LinearAlgebra._sqrt_real_2x2!)($y, $x);
  38.129 ns (0 allocations: 0 bytes)

julia> @btime $(LinearAlgebra.sqrt_quasitriu)($(F.T));
  17.354 μs (1 allocation: 3.25 KiB)

[musm]

This is a bugfix so this probably needs to be backported?

[sethaxen]

This is a bugfix so this probably needs to be backported?

No, this fixes bugs introduced in #39973. The new tests added here, which fail on master, pass on Julia 1.6.

[musm]

Merging, as this has been approved by @stevengj and has had enough time for any objections or comments since then.