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 😄

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.

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

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?

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

and now even works with (Furlong) units

So it is unit correct now.

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