Add diagview to obtain a view along a diagonal

[jishnub]

A function to obtain a view of a diagonal of a matrix is useful, and this is clearly being used widely within LinearAlgebra.

The implementation here iterates according to the IndexStyle of the array:

julia> using LinearAlgebra

julia> A = reshape(1:9, 3, 3)
3×3 reshape(::UnitRange{Int64}, 3, 3) with eltype Int64:
 1  4  7
 2  5  8
 3  6  9

julia> diagview(A,1)
2-element view(::UnitRange{Int64}, 4:4:8) with eltype Int64:
 4
 8

julia> T = Tridiagonal(1:3, 3:6, 4:6)
4×4 Tridiagonal{Int64, UnitRange{Int64}}:
 3  4  ⋅  ⋅
 1  4  5  ⋅
 ⋅  2  5  6
 ⋅  ⋅  3  6

julia> diagview(T,1)
3-element view(::Tridiagonal{Int64, UnitRange{Int64}}, StepRangeLen(CartesianIndex(1, 2), CartesianIndex(1, 1), 3)) with eltype Int64:
 4
 5
 6

Closes #30250

[mcabbott]

Seems OK, it's used a lot internally, and this is usually the only reason to want diagind.

However, it seems like the one-arg form should understand structured matrices, so that e.g. diagview(T) isa UnitRange. And likewise diagview(UpperTriangular(rand(4,4))) should probably act on the parent (which has linear indexing).

What diagview(UpperTriangular(rand(4,4)), 2) should do is less obvious... it could make an efficient view(::Matrix) in the useful case, but not for the -2 case, hence at the cost of type instability.

[jishnub]

Type-stability is indeed the reason I didn't specialize this for banded matrices. This parallels how diag doesn't copy the bands, but writes to a similar vector. While the single-argument method may be specialized, I wonder if this is particularly useful? Accessing any band other than the principal diagonal would require the 2-argument method. We would also want to ensure that diagview(A) and diagview(A,0) return the same type, so whether we specialize the one-argument method goes along with whether we want type-stability in the two-argument method.

[mcabbott]

We would also want to ensure that diagview(A) and diagview(A,0) return the same type

Why do we want this?

[jishnub]

Mainly for ease of use. Having code that works for diagview(A) but not for diagview(A,0) doesn't sound ideal.

[mcabbott]

I think you're suggesting someone may pass diagview(A) to a function which accepts only ::Vector, and may be surprised if they later change to diagview(A, 0). That seems unlikely to me, and sort-of backwards? If you ask for a view of some unknown AbstractMatrix, you generically expect some AbstractVector whose type is complicated... and if you accept that, then when someone passes in A::Diagonal, getting a simple Vector seems like it should be fine.

[jishnub]

I was thinking of use cases where diag may be swapped with diagview to make a call minimally allocating. E.g.:

julia> using LinearAlgebra

julia> A = rand(4,4)
4×4 Matrix{Float64}:
 0.805581  0.957599  0.252321   0.552123
 0.663896  0.821837  0.788979   0.110341
 0.654359  0.716732  0.0641013  0.0818827
 0.903501  0.174516  0.193872   0.933869

julia> @btime Bidiagonal(diag($A), diag($A,1), :U)
  82.477 ns (4 allocations: 176 bytes)
4×4 Bidiagonal{Float64, Vector{Float64}}:
 0.805581  0.957599   ⋅          ⋅ 
  ⋅        0.821837  0.788979    ⋅ 
  ⋅         ⋅        0.0641013  0.0818827
  ⋅         ⋅         ⋅         0.933869

julia> @btime Bidiagonal(diagview($A), diagview($A,1), :U)
  52.661 ns (2 allocations: 64 bytes)
4×4 Bidiagonal{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{StepRange{Int64, Int64}}, true}}:
 0.805581  0.957599   ⋅          ⋅ 
  ⋅        0.821837  0.788979    ⋅ 
  ⋅         ⋅        0.0641013  0.0818827
  ⋅         ⋅         ⋅         0.933869

The constructor only accepts AbstractVectors of the same type, so having diag work and diagview not work might be puzzling.

That said, having a direct view into the parent presents significant optimization opportunities, so I'm open to reconsidering the strict type-stability requirement.

[mcabbott]

Ok, fair point. With A::Matrix it will always be fine, but with UpperTriangular you would need to use diagview(U,0) to get two views of the same type:

julia> diagview(A::AbstractMatrix, k::Integer=0) = @view A[diagind(A, k, IndexStyle(A))]  # as in PR

julia> diagview(A::UpperTriangular) = diagview(parent(A))  # proposed optimisation

julia> U = UpperTriangular(A);

julia> Bidiagonal(diagview(U), diagview(U,1), :U)  # not the most informative error!
ERROR: MethodError: no method matching Bidiagonal(::SubArray{…}, ::SubArray{…}, ::Symbol)

julia> Bidiagonal(diagview(U,0), diagview(U,1), :U)  # like this is fine
4×4 Bidiagonal{Float64, SubArray{Float64, 1, UpperTriangular{Float64, Matrix{Float64}}, Tuple{StepRangeLen{CartesianIndex{2}, CartesianIndex{2}, CartesianIndex{2}, Int64}}, false}}:
 0.875899  0.396184   ⋅          ⋅ 
  ⋅        0.466412  0.175884    ⋅ 
  ⋅         ⋅        0.0872872  0.324497
  ⋅         ⋅         ⋅         0.142097

[jishnub]

I wonder if diagview should accept an argument to specify if the view should be translated to the parent? By default, it may return a view into the top-level array as in this PR, but one may opt into the alternate behavior of obtaining a band of the parent (or a view into one, e.g. for a BandedMatrix). This way, we retain the type-stability in the standard use case, but one doesn't need to sacrifice performance either.

[jishnub]

Gentle bump

[dkarrasch]

AFAIU, this is good for the time being, and (tempting) specializations could be added separately after careful discussion and consideration, right?

[longemen3000]

Just saw this PR, does it make sense to add diagview(x::Number) = x ?

[jishnub]

Is there motivation for this? To apply this function recursively, perhaps?