From db51022bce326833816b4505b17b15f9c82df3d2 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 14 Jul 2020 11:19:22 +0100 Subject: [PATCH 1/5] Stop transpose(A) \ lu(B)' from overwriting `A` MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This fixes the following bug I just found: ```julia julia> A = randn(5,5) 5×5 Array{Float64,2}: 0.485776 1.29655 0.0790172 0.66021 -1.49472 0.971676 -1.01382 0.630476 0.479027 -0.843428 -0.264609 0.392383 0.646729 0.510696 0.34913 -0.795944 -2.47709 -1.81052 -1.0947 -0.30381 -0.938873 2.16395 -1.33484 -0.745461 1.43709 julia> transpose(A) / lu(randn(5,5))' # should not change A 5×5 Adjoint{Float64,Array{Float64,2}}: 14.36 11.9797 -7.32563 1.87016 -16.2731 -18.4415 -13.7036 10.6455 -2.47396 19.9999 8.0401 8.31723 -5.16714 2.13261 -9.637 4.19849 3.9865 -1.98478 0.714778 -4.62445 -5.36059 -2.60991 0.917052 0.290281 4.62547 julia> A # but does! 5×5 Array{Float64,2}: 14.36 -18.4415 8.0401 4.19849 -5.36059 11.9797 -13.7036 8.31723 3.9865 -2.60991 -7.32563 10.6455 -5.16714 -1.98478 0.917052 1.87016 -2.47396 2.13261 0.714778 0.290281 -16.2731 19.9999 -9.637 -4.62445 4.62547 ``` Still need to add tests. --- stdlib/LinearAlgebra/src/lu.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 54d421441b9de..7f0345508ff96 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -432,11 +432,11 @@ end (/)(adjA::Adjoint{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) = adjoint(F.parent \ adjA.parent) function (/)(trA::Transpose{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(trA), eltype(F)) - return adjoint(ldiv!(F.parent, convert(AbstractVector{T}, conj(trA.parent)))) + return adjoint(ldiv!(F.parent, convert(AbstractVector{T}, conj!(copy(trA.parent))))) end function (/)(trA::Transpose{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(trA), eltype(F)) - return adjoint(ldiv!(F.parent, convert(AbstractMatrix{T}, conj(trA.parent)))) + return adjoint(ldiv!(F.parent, convert(AbstractMatrix{T}, conj!(copy(trA.parent))))) end function det(F::LU{T}) where T From 12e9b50a27395ee430f3632cd7f8bbd4782baa78 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 14 Jul 2020 11:25:41 +0100 Subject: [PATCH 2/5] Add tests for transpose(A) / lu(B)' --- stdlib/LinearAlgebra/test/lu.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index 1f5f2892947dc..048f2a7fe8005 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -360,4 +360,17 @@ end end end +@testset "transpose(A) / lu(B)' should not overwrite A (#36657)" begin + for elty in (Float16, Float64, ComplexF64) + A = randn(elty, 5, 5) + B = randn(elty, 5, 5) + C = copy(A) + a = randn(elty, 5) + c = copy(a) + @test transpose(A) / lu(B)' ≈ transpose(A) / B + @test transpose(a) / lu(B)' ≈ transpose(a) / B + @test A == C + @test a == c + end +end end # module TestLU From c9424fde957a9b4e46211b66c6bc462fa1cc7267 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 14 Jul 2020 11:31:14 +0100 Subject: [PATCH 3/5] Fix mixed types --- stdlib/LinearAlgebra/src/lu.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 7f0345508ff96..0d645f43b5c25 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -432,11 +432,11 @@ end (/)(adjA::Adjoint{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) = adjoint(F.parent \ adjA.parent) function (/)(trA::Transpose{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(trA), eltype(F)) - return adjoint(ldiv!(F.parent, convert(AbstractVector{T}, conj!(copy(trA.parent))))) + return adjoint(ldiv!(F.parent, conj!(AbstractVector{T}(trA.parent)))) end function (/)(trA::Transpose{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(trA), eltype(F)) - return adjoint(ldiv!(F.parent, convert(AbstractMatrix{T}, conj!(copy(trA.parent))))) + return adjoint(ldiv!(F.parent, conj!(AbstractMatrix{T}(trA.parent)))) end function det(F::LU{T}) where T From 7d2c4b35e2315ccc44f898a147279b21a0729a8b Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 14 Jul 2020 20:13:20 +0100 Subject: [PATCH 4/5] Update stdlib/LinearAlgebra/test/lu.jl Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/test/lu.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index 048f2a7fe8005..02ebb55b6e7f2 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -367,8 +367,8 @@ end C = copy(A) a = randn(elty, 5) c = copy(a) - @test transpose(A) / lu(B)' ≈ transpose(A) / B - @test transpose(a) / lu(B)' ≈ transpose(a) / B + @test transpose(A) / lu(B)' ≈ transpose(A) / B' + @test transpose(a) / lu(B)' ≈ transpose(a) / B' @test A == C @test a == c end From 5ade46f4f7d0773a9ca16b261358861ba8e5cca4 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 15 Jul 2020 19:36:21 +0100 Subject: [PATCH 5/5] Use copy_oftype --- stdlib/LinearAlgebra/src/lu.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 0d645f43b5c25..9262b517f5906 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -432,11 +432,11 @@ end (/)(adjA::Adjoint{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) = adjoint(F.parent \ adjA.parent) function (/)(trA::Transpose{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(trA), eltype(F)) - return adjoint(ldiv!(F.parent, conj!(AbstractVector{T}(trA.parent)))) + return adjoint(ldiv!(F.parent, conj!(copy_oftype(trA.parent, T)))) end function (/)(trA::Transpose{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(trA), eltype(F)) - return adjoint(ldiv!(F.parent, conj!(AbstractMatrix{T}(trA.parent)))) + return adjoint(ldiv!(F.parent, conj!(copy_oftype(trA.parent, T)))) end function det(F::LU{T}) where T