Fix triu/tril for partly initialized matrices #55312

jishnub · 2024-07-30T16:02:22Z

This fixes

julia> using LinearAlgebra, StaticArrays

julia> M = Matrix{BigInt}(undef, 2, 2); M[1,1] = M[2,2] = M[1,2] = 3;

julia> S = SizedMatrix{2,2}(M)
2×2 SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}} with indices SOneTo(2)×SOneTo(2):
   3     3
 #undef  3

julia> triu(S)
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex
   @ ./essentials.jl:907 [inlined]
 [2] getindex
   @ ~/.julia/packages/StaticArrays/MSJcA/src/SizedArray.jl:92 [inlined]
 [3] copyto_unaliased!
   @ ./abstractarray.jl:1086 [inlined]
 [4] copyto!(dest::SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}}, src::SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}})
   @ Base ./abstractarray.jl:1066
 [5] copymutable
   @ ./abstractarray.jl:1200 [inlined]
 [6] triu(M::SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/generic.jl:413
 [7] top-level scope
   @ REPL[11]:1

After this PR:

julia> triu(S)
2×2 SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}} with indices SOneTo(2)×SOneTo(2):
 3  3
 0  3

Only the indices that need to be copied are accessed, and the others are written to without being read.

~~Something similar needs to be done for the triu(A, k) and tril(A, k) methods as well.~~

dkarrasch · 2024-08-07T07:20:18Z

What about

triu([randn(2,2) for _ in 1:2, _ in 1:2])

? Matrices whose eltype doesn't admit a zero method, but whose elements admit one.

EDIT: Ouch, that seems to fail:

julia> triu([randn(2,2) for _ in 1:2, _ in 1:2])
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex
   @ ./essentials.jl:905 [inlined]
 [2] getindex
   @ ./array.jl:919 [inlined]
 [3] _zero
   @ ~/.julia/juliaup/julia-1.11.0-rc2+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:112 [inlined]
 [4] triu!(M::Matrix{Matrix{Float64}}, k::Int64)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.0-rc2+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:144
 [5] triu!
   @ ~/.julia/juliaup/julia-1.11.0-rc2+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:441 [inlined]
 [6] triu(M::Matrix{Matrix{Float64}})
   @ Main ./REPL[4]:3

As is introduces regression...

jishnub · 2024-08-07T09:03:39Z

Oops, triu! does require the matrix to be initialized. Need to figure out the order of triu! and copytrito!, largely because similar(::SymTridiagonal) returns a SymTridiagonal, but this doesn't support setindex!. Perhaps this should be modified to return a Tridiagonal instead?

dkarrasch · 2024-10-17T12:16:28Z

Is this ready?

jishnub · 2024-10-17T13:39:34Z

Yes, this is ready

Fixes https://github.com/JuliaLang/julia/issues/56134 After this, ```julia julia> using LinearAlgebra julia> A = hermitianpart(rand(4, 4)) 4×4 Hermitian{Float64, Matrix{Float64}}: 0.387617 0.277226 0.67629 0.60678 0.277226 0.894101 0.388416 0.489141 0.67629 0.388416 0.100907 0.619955 0.60678 0.489141 0.619955 0.452605 julia> B = UpperTriangular(A) 4×4 UpperTriangular{Float64, Hermitian{Float64, Matrix{Float64}}}: 0.387617 0.277226 0.67629 0.60678 ⋅ 0.894101 0.388416 0.489141 ⋅ ⋅ 0.100907 0.619955 ⋅ ⋅ ⋅ 0.452605 julia> B - B' 4×4 Matrix{Float64}: 0.0 0.277226 0.67629 0.60678 -0.277226 0.0 0.388416 0.489141 -0.67629 -0.388416 0.0 0.619955 -0.60678 -0.489141 -0.619955 0.0 ``` This preserves the band structure of the parent, if any: ```julia julia> U = UpperTriangular(Diagonal(ones(4))) 4×4 UpperTriangular{Float64, Diagonal{Float64, Vector{Float64}}}: 1.0 0.0 0.0 0.0 ⋅ 1.0 0.0 0.0 ⋅ ⋅ 1.0 0.0 ⋅ ⋅ ⋅ 1.0 julia> U - U' 4×4 Diagonal{Float64, Vector{Float64}}: 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ``` This doesn't fully work with partly initialized matrices, and would need #55312 for that. The abstract triangular methods now construct matrices using `similar(parent(U), size(U))` so that the destinations are fully mutable. ```julia julia> @invoke B::LinearAlgebra.AbstractTriangular - B'::LinearAlgebra.AbstractTriangular 4×4 Matrix{Float64}: 0.0 0.277226 0.67629 0.60678 -0.277226 0.0 0.388416 0.489141 -0.67629 -0.388416 0.0 0.619955 -0.60678 -0.489141 -0.619955 0.0 ``` --------- Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>

This fixes ```julia julia> using LinearAlgebra, StaticArrays julia> M = Matrix{BigInt}(undef, 2, 2); M[1,1] = M[2,2] = M[1,2] = 3; julia> S = SizedMatrix{2,2}(M) 2×2 SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}} with indices SOneTo(2)×SOneTo(2): 3 3 #undef 3 julia> triu(S) ERROR: UndefRefError: access to undefined reference Stacktrace: [1] getindex @ ./essentials.jl:907 [inlined] [2] getindex @ ~/.julia/packages/StaticArrays/MSJcA/src/SizedArray.jl:92 [inlined] [3] copyto_unaliased! @ ./abstractarray.jl:1086 [inlined] [4] copyto!(dest::SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}}, src::SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}}) @ Base ./abstractarray.jl:1066 [5] copymutable @ ./abstractarray.jl:1200 [inlined] [6] triu(M::SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}}) @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/generic.jl:413 [7] top-level scope @ REPL[11]:1 ``` After this PR: ```julia julia> triu(S) 2×2 SizedMatrix{2, 2, BigInt, 2, Matrix{BigInt}} with indices SOneTo(2)×SOneTo(2): 3 3 0 3 ``` Only the indices that need to be copied are accessed, and the others are written to without being read.

Fixes https://github.com/JuliaLang/julia/issues/56134 After this, ```julia julia> using LinearAlgebra julia> A = hermitianpart(rand(4, 4)) 4×4 Hermitian{Float64, Matrix{Float64}}: 0.387617 0.277226 0.67629 0.60678 0.277226 0.894101 0.388416 0.489141 0.67629 0.388416 0.100907 0.619955 0.60678 0.489141 0.619955 0.452605 julia> B = UpperTriangular(A) 4×4 UpperTriangular{Float64, Hermitian{Float64, Matrix{Float64}}}: 0.387617 0.277226 0.67629 0.60678 ⋅ 0.894101 0.388416 0.489141 ⋅ ⋅ 0.100907 0.619955 ⋅ ⋅ ⋅ 0.452605 julia> B - B' 4×4 Matrix{Float64}: 0.0 0.277226 0.67629 0.60678 -0.277226 0.0 0.388416 0.489141 -0.67629 -0.388416 0.0 0.619955 -0.60678 -0.489141 -0.619955 0.0 ``` This preserves the band structure of the parent, if any: ```julia julia> U = UpperTriangular(Diagonal(ones(4))) 4×4 UpperTriangular{Float64, Diagonal{Float64, Vector{Float64}}}: 1.0 0.0 0.0 0.0 ⋅ 1.0 0.0 0.0 ⋅ ⋅ 1.0 0.0 ⋅ ⋅ ⋅ 1.0 julia> U - U' 4×4 Diagonal{Float64, Vector{Float64}}: 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ``` This doesn't fully work with partly initialized matrices, and would need #55312 for that. The abstract triangular methods now construct matrices using `similar(parent(U), size(U))` so that the destinations are fully mutable. ```julia julia> @invoke B::LinearAlgebra.AbstractTriangular - B'::LinearAlgebra.AbstractTriangular 4×4 Matrix{Float64}: 0.0 0.277226 0.67629 0.60678 -0.277226 0.0 0.388416 0.489141 -0.67629 -0.388416 0.0 0.619955 -0.60678 -0.489141 -0.619955 0.0 ``` --------- Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>

Fixes https://github.com/JuliaLang/julia/issues/56134 After this, ```julia julia> using LinearAlgebra julia> A = hermitianpart(rand(4, 4)) 4×4 Hermitian{Float64, Matrix{Float64}}: 0.387617 0.277226 0.67629 0.60678 0.277226 0.894101 0.388416 0.489141 0.67629 0.388416 0.100907 0.619955 0.60678 0.489141 0.619955 0.452605 julia> B = UpperTriangular(A) 4×4 UpperTriangular{Float64, Hermitian{Float64, Matrix{Float64}}}: 0.387617 0.277226 0.67629 0.60678 ⋅ 0.894101 0.388416 0.489141 ⋅ ⋅ 0.100907 0.619955 ⋅ ⋅ ⋅ 0.452605 julia> B - B' 4×4 Matrix{Float64}: 0.0 0.277226 0.67629 0.60678 -0.277226 0.0 0.388416 0.489141 -0.67629 -0.388416 0.0 0.619955 -0.60678 -0.489141 -0.619955 0.0 ``` This preserves the band structure of the parent, if any: ```julia julia> U = UpperTriangular(Diagonal(ones(4))) 4×4 UpperTriangular{Float64, Diagonal{Float64, Vector{Float64}}}: 1.0 0.0 0.0 0.0 ⋅ 1.0 0.0 0.0 ⋅ ⋅ 1.0 0.0 ⋅ ⋅ ⋅ 1.0 julia> U - U' 4×4 Diagonal{Float64, Vector{Float64}}: 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 ``` This doesn't fully work with partly initialized matrices, and would need JuliaLang/julia#55312 for that. The abstract triangular methods now construct matrices using `similar(parent(U), size(U))` so that the destinations are fully mutable. ```julia julia> @invoke B::LinearAlgebra.AbstractTriangular - B'::LinearAlgebra.AbstractTriangular 4×4 Matrix{Float64}: 0.0 0.277226 0.67629 0.60678 -0.277226 0.0 0.388416 0.489141 -0.67629 -0.388416 0.0 0.619955 -0.60678 -0.489141 -0.619955 0.0 ``` --------- Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>

jishnub added linear algebra Linear algebra arrays [a, r, r, a, y, s] backport 1.11 Change should be backported to release-1.11 and removed backport 1.11 Change should be backported to release-1.11 labels Jul 30, 2024

jishnub requested a review from dkarrasch August 6, 2024 19:07

dkarrasch previously approved these changes Aug 7, 2024

View reviewed changes

jishnub marked this pull request as draft August 7, 2024 09:03

jishnub force-pushed the jishnub/generictriangularscaling branch from 51b27ef to 029ae61 Compare August 31, 2024 10:34

jishnub mentioned this pull request Oct 14, 2024

Specialize adding/subtracting mixed Upper/LowerTriangular #56149

Merged

jishnub added 5 commits October 14, 2024 11:55

Fix triu/tril for partly initialized matrices

c3e0c0f

Switch order of operations

39b5ecd

Adapt triu/tril with a band index

57a7ed0

Add similarmutable

2e4a69f

Merge triu/tril methods

7a3ef7e

jishnub force-pushed the jishnub/generictriangularscaling branch from 029ae61 to 7a3ef7e Compare October 14, 2024 07:31

jishnub added 2 commits October 14, 2024 20:12

Remove unused similarmutable

9daa8e9

Fast path for k==0

4329d07

jishnub marked this pull request as ready for review October 15, 2024 06:31

Merge branch 'master' into jishnub/generictriangularscaling

88006e3

dkarrasch added the merge me PR is reviewed. Merge when all tests are passing label Oct 17, 2024

dkarrasch approved these changes Oct 17, 2024

View reviewed changes

giordano merged commit f1990e2 into master Oct 17, 2024
8 checks passed

giordano deleted the jishnub/generictriangularscaling branch October 17, 2024 21:19

giordano removed the merge me PR is reviewed. Merge when all tests are passing label Oct 17, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix triu/tril for partly initialized matrices #55312

Fix triu/tril for partly initialized matrices #55312

jishnub commented Jul 30, 2024 •

edited

Loading

dkarrasch commented Aug 7, 2024 •

edited

Loading

jishnub commented Aug 7, 2024

dkarrasch commented Oct 17, 2024

jishnub commented Oct 17, 2024

Fix triu/tril for partly initialized matrices #55312

Fix triu/tril for partly initialized matrices #55312

Conversation

jishnub commented Jul 30, 2024 • edited Loading

dkarrasch commented Aug 7, 2024 • edited Loading

jishnub commented Aug 7, 2024

dkarrasch commented Oct 17, 2024

jishnub commented Oct 17, 2024

jishnub commented Jul 30, 2024 •

edited

Loading

dkarrasch commented Aug 7, 2024 •

edited

Loading