Adjoint for diagonal should be lazy

[jishnub]

Currently, adjoint for a Diagonal materializes the diagonal vector. This goes against the documentation of adjoint as a lazy wrapper.

This PR makes the result share memory with the original Diagonal.

julia> D = Diagonal(fill(2im, 2))
2×2 Diagonal{Complex{Int64}, Vector{Complex{Int64}}}:
 0+2im    ⋅  
   ⋅    0+2im

julia> D'
2×2 Diagonal{Complex{Int64}, Base.ReshapedArray{Complex{Int64}, 1, Adjoint{Complex{Int64}, Vector{Complex{Int64}}}, Tuple{}}}:
 0-2im    ⋅  
   ⋅    0-2im

The implementation isn't perfect, as

julia> D = Diagonal(fill(2im, 2))
2×2 Diagonal{Complex{Int64}, Vector{Complex{Int64}}}:
 0+2im    ⋅  
   ⋅    0+2im

julia> D'
2×2 Diagonal{Complex{Int64}, Base.ReshapedArray{Complex{Int64}, 1, Adjoint{Complex{Int64}, Vector{Complex{Int64}}}, Tuple{}}}:
 0-2im    ⋅  
   ⋅    0-2im

julia> (D')'
2×2 Diagonal{Complex{Int64}, Base.ReshapedArray{Complex{Int64}, 1, Adjoint{Complex{Int64}, Base.ReshapedArray{Complex{Int64}, 1, Adjoint{Complex{Int64}, Vector{Complex{Int64}}}, Tuple{}}}, Tuple{}}}:
 0+2im    ⋅  
   ⋅    0+2im

and, ideally, (D')' would just return D. However, at least this makes the adjoint operation lazy and allocation-free.

[dkarrasch]

SGTM. Does that affect multiplication/solving performance somehow?

[dkarrasch]

The other option would be, of course, to make D' = Adjoint(D), but then we would need to add a couple of multiplication and solve methods (which I removed at some point because one would not get an Adjoint via adjoint, irrespective of the eltype).

[ViralBShah]

Good to merge?

[jishnub]

Not yet, I'll check the performance implication

[jishnub]

Performance for some simple operations seems almost unaffected. If anything, solve seems slightly improved:
On master

julia> D = Diagonal([1:1000;]*im)';

julia> v = rand(ComplexF64, 1000);

julia> @btime $D * $D;
  1.611 μs (1 allocation: 15.75 KiB)

julia> @btime $D * $v;
  1.672 μs (1 allocation: 15.75 KiB)

julia> @btime $D \ $v;
  31.136 μs (1 allocation: 15.75 KiB)

julia> @btime $D \ $D;
  23.815 μs (1 allocation: 15.75 KiB)

julia> @btime $D \ $(Matrix(D));
  22.762 ms (2 allocations: 15.26 MiB)

This PR

julia> @btime $D * $D;
  1.610 μs (1 allocation: 15.75 KiB)

julia> @btime $D * $v;
  1.766 μs (1 allocation: 15.75 KiB)

julia> @btime $D \ $v;
  25.867 μs (1 allocation: 15.75 KiB)

julia> @btime $D \ $D;
  7.162 μs (1 allocation: 15.75 KiB)

julia> @btime $D \ $(Matrix(D));
  22.506 ms (2 allocations: 15.26 MiB)

[dkarrasch]

For the double adjoint, couldn't you add a method for ::Diagonal{<:Any,<:ReshapedArray{<:Adjoint}} and make it do the ideal thing? And mitigate the issue further via

adjoint(D::Diagonal{<:Real}) = transpose(D)

? Then there should remain a smaller niche for this "ugly" behavior.

[jishnub]

Updated, now

julia> D = Diagonal([1,2,3]*im)
3×3 Diagonal{Complex{Int64}, Vector{Complex{Int64}}}:
 0+1im    ⋅      ⋅  
   ⋅    0+2im    ⋅  
   ⋅      ⋅    0+3im

julia> (D')'
3×3 Diagonal{Complex{Int64}, Vector{Complex{Int64}}}:
 0+1im    ⋅      ⋅  
   ⋅    0+2im    ⋅  
   ⋅      ⋅    0+3im

The real case already does the right thing, so no special-casing is necessary.