Improve performance of svd and eigen of Diagonals

[dkarrasch]

This reduces lots of intermediate array allocations, and now even works with (Furlong) units. I'll add a test for that case later.

[dkarrasch]

This is a wonderful example why splitting LinearAlgebra and SparseArrays was really useful. See the benchmarks:

julia> D = Diagonal(randn(1024));

julia> @btime svd($D);
  16.898 ms (10275 allocations: 866.75 KiB) # old behavior with sparse U and Vt

julia> @btime svd($D);
  11.371 ms (2076 allocations: 40.38 MiB) # before this PR: dense output + temporary allocations

julia> @btime svd($D);
  1.506 ms (8 allocations: 16.02 MiB) # this PR: dense output + seemingly no unnecessary allocations

While sparse output is good from a memory perspective, writing into an initially empty sparse matrix is slow. I think the same is true for the multiplications listed in 14154fc#commitcomment-64103782. Honestly, I hadn't seen that coming. 😅 So, I believe the "new" defaults (similar(::StructuredMatrix, T, size) = zeros() instead of spzeros) are really not bad, plus the old behavior is easy to recover.

Pinging some involved people to make sure they join me in my little github celebration 🎉 : @vtjnash @KristofferC @ViralBShah 😄

[dkarrasch]

This also fixes an issue in sorted eigen of Diagonals, where the sort! was acting on D.diag directly. I'll prepare a minimal fix suitable for backporting.

[KristofferC]

I tried to backport this to 1.6 but got a test error:

julia> E = eigen(Du)
Eigen{Furlong{0, Float64}, Furlong{1, Float64}, Matrix{Furlong{0, Float64}}, Vector{Furlong{1, Float64}}}
values:
3-element Vector{Furlong{1, Float64}}:
  Furlong{1, Float64}(-0.5153917252601611)
 Furlong{1, Float64}(-0.23809233488008882)
   Furlong{1, Float64}(0.1571713328338422)
vectors:
3×3 Matrix{Furlong{0, Float64}}:
 Furlong{0, Float64}(1.0)  Furlong{0, Float64}(0.0)  Furlong{0, Float64}(0.0)
 Furlong{0, Float64}(0.0)  Furlong{0, Float64}(1.0)  Furlong{0, Float64}(0.0)
 Furlong{0, Float64}(0.0)  Furlong{0, Float64}(0.0)  Furlong{0, Float64}(1.0)

julia> @test Matrix(E) == Du
Error During Test at REPL[55]:1
  Test threw exception
  Expression: Matrix(E) == Du
  AssertionError: p == q
  Stacktrace:
    [1] Furlong{1, Float64}(x::Furlong{0, Float64})
      @ Main.Furlongs ~/julia/test/testhelpers/Furlongs.jl:20
    [2] convert(#unused#::Type{Furlong{1, Float64}}, x::Furlong{0, Float64})
...

at

julia/stdlib/LinearAlgebra/test/diagonal.jl

Line 427 in ac1d693

@test map(x -> x.val, Matrix(F)) ≈ map(x -> x.val, Du)

[dkarrasch]

I can reproduce. I think the diagonal multiplication code must have changed between versions, and on v1.6 it cannot multiply vecs::Matrix{Furlong{0}} * Diagonal(vals::Vector{Furlong{1}}. It tries to copy vecs into a matrix of the result type, which impossible. Long story short, it fails for reasons outside of this commit, so can we still backport but omit the unitful tests?

[oscardssmith]

why are we back porting? isn't this just a performance improvement?

[KristofferC]

and now even works with (Furlong) units

So it is unit correct now.

[dkarrasch]

First of all, it contains a fix for the aliasing of D.diag in sorted eigendecompositions.