Skip to content

Commit cc74d24

Browse files
authored
Inplace transpose for unit Triangular may skip diagonal (JuliaLang#53101)
Since the diagonal elements of a `UnitUpperTriangular` are given by `onelement`, these should be unchanged under `transpose/adjoint`, and we don't need to access these elements in the parent array when performing in-place operations. Fixes ```julia julia> using LinearAlgebra julia> M = Matrix{BigFloat}(undef, 2, 2); julia> M[1,2] = 3; julia> U = UnitUpperTriangular(M) 2×2 UnitUpperTriangular{BigFloat, Matrix{BigFloat}}: 1.0 3.0 ⋅ 1.0 julia> transpose!(U) ERROR: UndefRefError: access to undefined reference Stacktrace: [1] getindex @ ./essentials.jl:882 [inlined] [2] getindex @ ./array.jl:915 [inlined] [3] copytri! @ ~/packages/julias/julia-latest/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:414 [inlined] [4] transpose!(A::UnitUpperTriangular{BigFloat, Matrix{BigFloat}}) @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:470 [5] top-level scope @ REPL[5]:1 ``` After this PR: ```julia julia> transpose!(U) 2×2 UnitLowerTriangular{BigFloat, Matrix{BigFloat}}: 1.0 ⋅ 3.0 1.0 ```
1 parent cb9b00d commit cc74d24

File tree

2 files changed

+27
-7
lines changed

2 files changed

+27
-7
lines changed

stdlib/LinearAlgebra/src/triangular.jl

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -465,13 +465,13 @@ transpose(A::UnitLowerTriangular) = UnitUpperTriangular(transpose(A.data))
465465
transpose(A::UnitUpperTriangular) = UnitLowerTriangular(transpose(A.data))
466466

467467
transpose!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L', false, true))
468-
transpose!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L', false, true))
468+
transpose!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L', false, false))
469469
transpose!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U', false, true))
470-
transpose!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U', false, true))
470+
transpose!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U', false, false))
471471
adjoint!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L' , true, true))
472-
adjoint!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L' , true, true))
472+
adjoint!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L' , true, false))
473473
adjoint!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U' , true, true))
474-
adjoint!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U' , true, true))
474+
adjoint!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U' , true, false))
475475

476476
diag(A::LowerTriangular) = diag(A.data)
477477
diag(A::UnitLowerTriangular) = fill(oneunit(eltype(A)), size(A,1))

stdlib/LinearAlgebra/test/triangular.jl

Lines changed: 23 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ debug && println("Test basic type functionality")
2626
@test LowerTriangular(randn(3, 3)) |> t -> [size(t, i) for i = 1:3] == [size(Matrix(t), i) for i = 1:3]
2727

2828
# The following test block tries to call all methods in base/linalg/triangular.jl in order for a combination of input element types. Keep the ordering when adding code.
29-
for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat}, Int)
29+
@testset for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat}, Int)
3030
# Begin loop for first Triangular matrix
3131
for (t1, uplo1) in ((UpperTriangular, :U),
3232
(UnitUpperTriangular, :U),
@@ -998,9 +998,12 @@ end
998998

999999
@testset "arithmetic with partly uninitialized matrices" begin
10001000
@testset "$(typeof(A))" for A in (Matrix{BigFloat}(undef,2,2), Matrix{Complex{BigFloat}}(undef,2,2)')
1001-
A[1,1] = A[2,2] = A[2,1] = 4
1001+
A[2,1] = eltype(A) <: Complex ? 4 + 3im : 4
10021002
B = Matrix{eltype(A)}(undef, size(A))
10031003
for MT in (LowerTriangular, UnitLowerTriangular)
1004+
if MT == LowerTriangular
1005+
A[1,1] = A[2,2] = eltype(A) <: Complex ? 4 + 3im : 4
1006+
end
10041007
L = MT(A)
10051008
B .= 0
10061009
copyto!(B, L)
@@ -1009,13 +1012,23 @@ end
10091012
@test 2\L == 2\B
10101013
@test real(L) == real(B)
10111014
@test imag(L) == imag(B)
1015+
if A isa Matrix
1016+
@test transpose!(MT(copy(A))) == transpose(L)
1017+
@test adjoint!(MT(copy(A))) == adjoint(L)
1018+
else
1019+
@test_broken transpose!(MT(copy(A))) == transpose(L)
1020+
@test_broken adjoint!(MT(copy(A))) == adjoint(L)
1021+
end
10121022
end
10131023
end
10141024

10151025
@testset "$(typeof(A))" for A in (Matrix{BigFloat}(undef,2,2), Matrix{Complex{BigFloat}}(undef,2,2)')
1016-
A[1,1] = A[2,2] = A[1,2] = 4
1026+
A[1,2] = eltype(A) <: Complex ? 4 + 3im : 4
10171027
B = Matrix{eltype(A)}(undef, size(A))
10181028
for MT in (UpperTriangular, UnitUpperTriangular)
1029+
if MT == UpperTriangular
1030+
A[1,1] = A[2,2] = eltype(A) <: Complex ? 4 + 3im : 4
1031+
end
10191032
U = MT(A)
10201033
B .= 0
10211034
copyto!(B, U)
@@ -1024,6 +1037,13 @@ end
10241037
@test 2\U == 2\B
10251038
@test real(U) == real(B)
10261039
@test imag(U) == imag(B)
1040+
if A isa Matrix
1041+
@test transpose!(MT(copy(A))) == transpose(U)
1042+
@test adjoint!(MT(copy(A))) == adjoint(U)
1043+
else
1044+
@test_broken transpose!(MT(copy(A))) == transpose(U)
1045+
@test_broken adjoint!(MT(copy(A))) == adjoint(U)
1046+
end
10271047
end
10281048
end
10291049
end

0 commit comments

Comments
 (0)