Change diag(A::StructuredMatrix[, k=0]) #24324
Conversation
base/linalg/bidiag.jl
Outdated
| return M.uplo == 'U' ? M.ev : zeros(T, size(M,1)-1) | ||
| elseif n == -1 | ||
| return M.uplo == 'L' ? M.ev : zeros(T, size(M,1)-1) | ||
| return copy(M.dv) |
There was a problem hiding this comment.
Should change this to copy!(similar(M.dv, length(M.dv)), M.dv) to make certain we return the same type, since we then call the same similar method no matter what n is.
fredrikekre
force-pushed
the
fe/diagfixes
branch
from
8c47b5c to
b7c13a5
Compare
| return M.uplo == 'U' ? M.ev : zeros(T, size(M,1)-1) | ||
| elseif n == -1 | ||
| return M.uplo == 'L' ? M.ev : zeros(T, size(M,1)-1) | ||
| return copy!(similar(M.dv, length(M.dv)), M.dv) |
There was a problem hiding this comment.
Could the length(M.dv) be omitted? (Edit: Guarding against the length specification modifying the return type?)
There was a problem hiding this comment.
Moreover, why not simply copy(... rather than copy!(similar(M.dv, length(M.dv)), M.dv)? Does the return type otherwise depend upon k in some cases (e.g. for ranges)?
There was a problem hiding this comment.
Guarding against the length specification modifying the return type?
Yes, that was the idea, to guarantee that we take the same path no matter which diagonal was requested.
Moreover, why not simply copy(... rather than copy!(similar(M.dv, length(M.dv)), M.dv)? Does the return type otherwise depend upon k in some cases (e.g. for ranges)?
Same comment applies; guarantee the same code path. Also copy(r::AbstractRange) = r so that will not yield a Vector{Int} as similar does.
julia> copy(1:4)
1:4I don't think the copy!(similar(...), ...) uses are that terrible -- it is pretty clear what is going on IMO.
There was a problem hiding this comment.
Sounds good then! I might consider adding comments to make the intention clear for future readers; to someone not acquainted with the subtleties involved, the intention may be far from clear :).
Sacha0
left a comment
Sacha0
left a comment
There was a problem hiding this comment.
If my guesses above at the reasons for the particular copy!(similar(... invocations are correct, then lgtm! :)
- diag(A::StructuredMatrix[, k=0]) now return a new vector,
such that it does not alias the wrapped vector. This now
behaves the same for eg Matrix, where we get a new vector.
It is also consistent with getindex (diag is essentially a
special case of getindex).
- diag(A::StructuredMatrix{T,VectorType}[, k=0]) now call
similar(..., ::Int) on every branch to make sure the same
vector type is always returned, this is usually of type
VectorType.
fredrikekre
force-pushed
the
fe/diagfixes
branch
from
b7c13a5 to
393e88f
Compare
| function diag(M::Bidiagonal{T}, n::Integer=0) where T | ||
| function diag(M::Bidiagonal, n::Integer=0) | ||
| # every branch call similar(..., ::Int) to make sure the | ||
| # same vector type is returned independent of n |
My take on JuliaLang/LinearAlgebra.jl#478
Changes:
diag(A::StructuredMatrix[, k=0])now return a new vector, such that it does not alias the wrapped vector. This now behaves the same for egMatrix, where we get a new vector. It is also consistent withgetindex(diagis essentially a special case ofgetindex).diag(A::StructuredMatrix{T,VectorType}[, k=0])always return a vector of typeVectorType, i.e. the same type which the diagonals are stored as (actually we returnsimilar(::VectorType, k::Int), but that is usually of typeVectorType).This is problematic for wrappedEDIT: The latest revision does actually work forUnitRanges but that is already pretty broken for this case.UnitRanges also, sincesimilar(1:4)return anVector{Int}🎉fix JuliaLang/LinearAlgebra.jl#478