Skip to content
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.

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

Access only stored inds in copy for strided AbstractTriangular #52907

Merged
merged 2 commits into from
Jan 31, 2024

Conversation

jishnub
Copy link
Contributor

@jishnub jishnub commented Jan 15, 2024

This PR fixes the following:
on master

julia> M = Matrix{Complex{BigFloat}}(undef, 2, 2);

julia> M[1,1] = M[2,2] = M[1,2] = 2;

julia> U = UpperTriangular(view(M, :, :)) # need a wrapper, as the error doesn't happen for `Array`s
2×2 UpperTriangular{Complex{BigFloat}, SubArray{Complex{BigFloat}, 2, Matrix{Complex{BigFloat}}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}}, true}}:
 2.0+0.0im  2.0+0.0im
           2.0+0.0im

julia> copy(U)
ERROR: UndefRefError: access to undefined reference
Stacktrace:
  [1] getindex
    @ ./essentials.jl:817 [inlined]
  [2] getindex
    @ ./array.jl:912 [inlined]
  [3] macro expansion
    @ ./multidimensional.jl:919 [inlined]
  [4] macro expansion
    @ ./cartesian.jl:64 [inlined]
  [5] _unsafe_getindex!
    @ ./multidimensional.jl:914 [inlined]
  [6] _unsafe_getindex(::IndexLinear, ::Matrix{Complex{BigFloat}}, ::Base.Slice{Base.OneTo{Int64}}, ::Base.Slice{Base.OneTo{Int64}})
    @ Base ./multidimensional.jl:905
  [7] _getindex
    @ Base ./multidimensional.jl:891 [inlined]
  [8] getindex
    @ Base ./abstractarray.jl:1307 [inlined]
  [9] copy
    @ Base ./subarray.jl:73 [inlined]
 [10] copy(A::UpperTriangular{Complex{BigFloat}, SubArray{Complex{BigFloat}, 2, Matrix{Complex{}}, Tuple{Base.Slice{}, Base.Slice{}}, true}})
    @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:47
 [11] top-level scope
    @ REPL[36]:1
Some type information was truncated. Use `show(err)` to see complete types

This PR:

julia> copy(U) # doesn't access the indices corresponding to structural zeros
2×2 UpperTriangular{Complex{BigFloat}, Matrix{Complex{BigFloat}}}:
 2.0+0.0im  2.0+0.0im
           2.0+0.0im

This also improves performance in isbits cases:

julia> M = Matrix{ComplexF64}(undef, 1000, 1000);

julia> L = LowerTriangular(M);

julia> @btime copy($L);
  1.234 ms (3 allocations: 15.26 MiB) # master
  843.973 μs (3 allocations: 15.26 MiB) # PR

@jishnub jishnub added linear algebra Linear algebra arrays [a, r, r, a, y, s] labels Jan 15, 2024
@jishnub
Copy link
Contributor Author

jishnub commented Jan 15, 2024

From the test failures, it looks like this would need #52528

@jishnub jishnub requested a review from dkarrasch January 26, 2024 05:51
@jishnub
Copy link
Contributor Author

jishnub commented Jan 28, 2024

Gentle bump

@dkarrasch
Copy link
Member

Seems like there's an issue with transposed triangular matrices?

@jishnub jishnub force-pushed the jishnub/copytristrided branch from 9946353 to 4bc03e4 Compare January 29, 2024 05:53
@jishnub
Copy link
Contributor Author

jishnub commented Jan 29, 2024

Could you provide an example? The following appears to work:

julia> A = Matrix{Complex{BigFloat}}(undef,2,2);

julia> A[1,1] = A[2,2] = A[1,2] = 4;

julia> L = LowerTriangular(transpose(A))
2×2 LowerTriangular{Complex{BigFloat}, Transpose{Complex{BigFloat}, Matrix{Complex{BigFloat}}}}:
 4.0+0.0im          
 4.0+0.0im  4.0+0.0im

julia> copy(L)
2×2 LowerTriangular{Complex{BigFloat}, Matrix{Complex{BigFloat}}}:
 4.0+0.0im          
 4.0+0.0im  4.0+0.0im

There are indeed cases where the fallback copy implementation is hit, which has issues with partly initialized matrices. However, the fix for that should probably go to copy for the parent types, such as transpose.

@dkarrasch
Copy link
Member

I mean, the tests were failing in triangular.jl on all platforms.

@jishnub
Copy link
Contributor Author

jishnub commented Jan 29, 2024

Oh, I had missed that somehow. The specific test that's failing should be fixed by #53101

@dkarrasch dkarrasch changed the title Access only stored inds in copy for strided AbstractTriangular Access only stored inds in copy for strided AbstractTriangular Jan 31, 2024
@dkarrasch dkarrasch added the merge me PR is reviewed. Merge when all tests are passing label Jan 31, 2024
@jishnub jishnub force-pushed the jishnub/copytristrided branch from a32190e to 1731cb5 Compare January 31, 2024 14:24
@IanButterworth IanButterworth merged commit 08e3c2e into master Jan 31, 2024
7 checks passed
@IanButterworth IanButterworth deleted the jishnub/copytristrided branch January 31, 2024 20:38
@IanButterworth IanButterworth removed the merge me PR is reviewed. Merge when all tests are passing label Jan 31, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
arrays [a, r, r, a, y, s] linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants