Fixed segfault when adding OffsetArray and UniformScaling #38544
Conversation
| @@ -215,17 +215,17 @@ end | |||
| function (+)(A::AbstractMatrix, J::UniformScaling) | |||
| checksquare(A) | |||
| B = copy_oftype(A, Base._return_type(+, Tuple{eltype(A), typeof(J)})) | |||
| @inbounds for i in axes(A, 1) | |||
| B[i,i] += J | |||
| for i in max(first(axes(A, 1)), first(axes(A, 2))):min(last(axes(A, 1)), last(axes(A, 2))) | |||
There was a problem hiding this comment.
Could this just be for i in intersect(axes(A)...)? Or less cryptically intersect(axes(A,1), axes(A,2)).
There was a problem hiding this comment.
Yes, intersect is exactly the needed function:
Line 824 in 809f27c
By the way, we can remove checksquare(A) test now. Is a new pull request needed to examine this proposal?
There was a problem hiding this comment.
For the intersect thing just update the PR. Removing checksquare seems like it would merit discussion in a separate PR.
There was a problem hiding this comment.
Although I guess this PR does raise the question whether this is really the right thing to do. So far, we've basically been looking at UniformScaling as an automatically resizing matrix, so it's not entirely clear what it should do with an OffsetArray. @timholy thoughts?
There was a problem hiding this comment.
@fp4code could you elaborate on the use case where you wanted this operation? That could perhaps guide intuition as to what it should do.
There was a problem hiding this comment.
There is two choices.
- extend +I as "fill the diagonal i==j with ones" ; then the constraint "checksquare" can be removed
- restrict +I to square matrices algebra, and indeed, checksquare needs to be rewritten in order to control the identity of axes
There was a problem hiding this comment.
For the intersect thing just update the PR...
Please can you help me about how to update a PR here? With just a new commit, or do you want a rebase?
There was a problem hiding this comment.
Worth noting that A = Matrix(I, 3, 4) does work. But it is a bit more strange for A + I to work, because normally you'd expect this to make sense:
v = rand(4)
(A + I)*v == A*v + v # RHS is an error
This property would also fail for an A with mismatched offsets. I guess that argues for option 2, making checksquare(A) check the axes. In which case there is no need to alter for i in axes(A, 1).
Is checksquare used in other contexts where allowing different index ranges on left & right would make sense?
There was a problem hiding this comment.
Please can you help me about how to update a PR here? With just a new commit, or do you want a rebase?
It doesn't really matter since the PR is likely to get squashed upon merge. Might be easiest to just add a commit.
There was a problem hiding this comment.
Is
checksquareused in other contexts where allowing different index ranges on left & right would make sense?
Not really. But checksquare is problematic now, because it is used to return the size and not the common axis range. Worst, the multiple argument checksquare, used in lapack.jl, is not efficient. See here for example
julia/stdlib/LinearAlgebra/src/lapack.jl
Line 6152 in 3074306
The multiple arguments
checksquare
is called, creating a vector by several inefficients push! calls. After this call, lapack function is checking that all matrices have the same dimensions. Well, the checksquare(A...) need to be replaced in lapack.jl by a more efficient function, get_common_size_of_square_matrices(A...). And here in uniformscaling.jl checksquare(A) could be replaced by r = get_range_of_square_matrix(A), replacing for i in axes(A, 1) by for i in r.
Fixes #38543.