Skip to content

Fixed segfault when adding OffsetArray and UniformScaling #38544

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Apr 3, 2021
Merged

Fixed segfault when adding OffsetArray and UniformScaling #38544

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a84482d
Fixed segfault when adding OffsetArray and UniformScaling
fp4code Nov 23, 2020
58b1c35
Using intersect function (suggested by @mcabbott)
fp4code Nov 30, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions stdlib/LinearAlgebra/src/uniformscaling.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)))
Copy link
Contributor

@mcabbott mcabbott Nov 23, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could this just be for i in intersect(axes(A)...)? Or less cryptically intersect(axes(A,1), axes(A,2)).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, intersect is exactly the needed function:

julia/base/range.jl

Line 824 in 809f27c

intersect(r::AbstractUnitRange{<:Integer}, s::AbstractUnitRange{<:Integer}) = max(first(r),first(s)):min(last(r),last(s))

By the way, we can remove checksquare(A) test now. Is a new pull request needed to examine this proposal?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the intersect thing just update the PR. Removing checksquare seems like it would merit discussion in a separate PR.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@fp4code could you elaborate on the use case where you wanted this operation? That could perhaps guide intuition as to what it should do.

Copy link
Contributor Author

@fp4code fp4code Nov 30, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is two choices.

  1. extend +I as "fill the diagonal i==j with ones" ; then the constraint "checksquare" can be removed
  2. restrict +I to square matrices algebra, and indeed, checksquare needs to be rewritten in order to control the identity of axes

Copy link
Contributor Author

@fp4code fp4code Nov 30, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Contributor

@mcabbott mcabbott Nov 30, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

@fp4code fp4code Nov 30, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is checksquare used 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

n, nt, nq, nz = checksquare(S, T, Q, Z)

The multiple arguments checksquare

julia/stdlib/LinearAlgebra/src/LinearAlgebra.jl

Line 227 in 3074306

function checksquare(A...)
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.

@inbounds B[i,i] += J
end
return B
end

function (-)(J::UniformScaling, A::AbstractMatrix)
checksquare(A)
B = convert(AbstractMatrix{Base._return_type(+, Tuple{eltype(A), typeof(J)})}, -A)
@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)))
@inbounds B[i,i] += J
end
return B
end
Expand Down
10 changes: 10 additions & 0 deletions stdlib/LinearAlgebra/test/uniformscaling.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@ using Test, LinearAlgebra, Random, SparseArrays
const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :Quaternions) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Quaternions.jl"))
using .Main.Quaternions
isdefined(Main, :OffsetArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl"))
using .Main.OffsetArrays

Random.seed!(123)

Expand Down Expand Up @@ -504,4 +506,12 @@ end
end
end

@testset "offset arrays" begin
A = OffsetArray(zeros(4,4), -1:2, 0:3)
@test sum(I + A) ≈ 3.0
@test sum(A + I) ≈ 3.0
@test sum(I - A) ≈ 3.0
@test sum(A - I) ≈ -3.0
end

end # module TestUniformscaling