From 6f33e52722f17d555881c7b768c31a418e220532 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 25 Jul 2023 20:41:12 +0200 Subject: [PATCH 1/4] Fix rdiv of complex lhs by real factorizations --- stdlib/LinearAlgebra/src/factorization.jl | 14 ++++++++------ stdlib/LinearAlgebra/test/hessenberg.jl | 3 ++- stdlib/LinearAlgebra/test/lu.jl | 9 +++++++++ 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index 8c35a23e6b6d5..7c9ff5dbdc313 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -151,17 +151,19 @@ end (\)(F::AdjointFactorization{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = @invoke \(F::typeof(F), B::VecOrMat) -function (/)(B::VecOrMat{Complex{T}}, F::Factorization{T}) where {T<:BlasReal} +function (/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}}}, F::Factorization{T}) where {T<:BlasReal} require_one_based_indexing(B) x = rdiv!(copy(reinterpret(T, B)), F) return copy(reinterpret(Complex{T}, x)) end # don't do the reinterpretation for [Adjoint/Transpose]Factorization -(/)(B::VecOrMat{Complex{T}}, F::TransposeFactorization{T}) where {T<:BlasReal} = - conj!(adjoint(parent(F)) \ conj.(B)) -(/)(B::VecOrMat{Complex{T}}, F::AdjointFactorization{T}) where {T<:BlasReal} = - @invoke /(B::VecOrMat{Complex{T}}, F::Factorization{T}) - +(/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}}}, F::TransposeFactorization{T}) where {T<:BlasReal} = + @invoke /(B::AbstractMatrix, F::Factorization) +(/)(B::Matrix{Complex{T}}, F::AdjointFactorization{T}) where {T<:BlasReal} = + @invoke /(B::AbstractMatrix, F::Factorization) +(/)(B::Adjoint{Complex{T},Vector{Complex{T}}}, F::AdjointFactorization{T}) where {T<:BlasReal} = + (F' \ B')' + function (\)(F::Factorization, B::AbstractVecOrMat) require_one_based_indexing(B) TFB = typeof(oneunit(eltype(F)) \ oneunit(eltype(B))) diff --git a/stdlib/LinearAlgebra/test/hessenberg.jl b/stdlib/LinearAlgebra/test/hessenberg.jl index 61e498211ca7b..62d10487b9b45 100644 --- a/stdlib/LinearAlgebra/test/hessenberg.jl +++ b/stdlib/LinearAlgebra/test/hessenberg.jl @@ -178,8 +178,9 @@ let n = 10 @test H \ B ≈ A \ B ≈ H \ complex(B) @test (H - I) \ B ≈ (A - I) \ B @test (H - (3+4im)I) \ B ≈ (A - (3+4im)I) \ B - @test b' / H ≈ b' / A ≈ complex.(b') / H + @test b' / H ≈ b' / A ≈ complex(b') / H @test B' / H ≈ B' / A ≈ complex(B') / H + @test b' / H' ≈ complex(b)' / H' @test B' / (H - I) ≈ B' / (A - I) @test B' / (H - (3+4im)I) ≈ B' / (A - (3+4im)I) @test (H - (3+4im)I)' \ B ≈ (A - (3+4im)I)' \ B diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index 25a75e13233f0..744cc543b993e 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -391,6 +391,15 @@ end B = randn(elty, 5, 5) @test rdiv!(transform(A), transform(lu(B))) ≈ transform(C) / transform(B) end + for elty in (Float32, Float64, ComplexF64), transF in (identity, transpose), + transB in (transpose, adjoint), transT in (identity, complex) + A = randn(elty, 5, 5) + F = lu(A) + b = randn(transT(elty), 5) + @test rdiv!(transB(copy(b)), transF(F)) ≈ transB(b) / transF(F) ≈ transB(b) / transF(A) + B = randn(transT(elty), 5, 5) + @test rdiv!(copy(B), transF(F)) ≈ B / transF(F) ≈ B / transF(A) + end end @testset "transpose(A) / lu(B)' should not overwrite A (#36657)" begin From d9d4e1f6593bf9fb0d0f0edb5f2df123010fd35c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 25 Jul 2023 20:43:43 +0200 Subject: [PATCH 2/4] move code --- stdlib/LinearAlgebra/src/factorization.jl | 56 +++++++++++------------ stdlib/LinearAlgebra/src/hessenberg.jl | 1 - stdlib/LinearAlgebra/src/lu.jl | 2 - 3 files changed, 28 insertions(+), 31 deletions(-) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index 7c9ff5dbdc313..bd93fccadb0ab 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -137,6 +137,12 @@ function Base.show(io::IO, ::MIME"text/plain", x::TransposeFactorization) show(io, MIME"text/plain"(), parent(x)) end +function (\)(F::Factorization, B::AbstractVecOrMat) + require_one_based_indexing(B) + TFB = typeof(oneunit(eltype(F)) \ oneunit(eltype(B))) + ldiv!(F, copy_similar(B, TFB)) +end +(\)(F::TransposeFactorization, B::AbstractVecOrMat) = conj!(adjoint(F.parent) \ conj.(B)) # With a real lhs and complex rhs with the same precision, we can reinterpret # the complex rhs as a real rhs with twice the number of columns or rows function (\)(F::Factorization{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} @@ -151,34 +157,6 @@ end (\)(F::AdjointFactorization{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = @invoke \(F::typeof(F), B::VecOrMat) -function (/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}}}, F::Factorization{T}) where {T<:BlasReal} - require_one_based_indexing(B) - x = rdiv!(copy(reinterpret(T, B)), F) - return copy(reinterpret(Complex{T}, x)) -end -# don't do the reinterpretation for [Adjoint/Transpose]Factorization -(/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}}}, F::TransposeFactorization{T}) where {T<:BlasReal} = - @invoke /(B::AbstractMatrix, F::Factorization) -(/)(B::Matrix{Complex{T}}, F::AdjointFactorization{T}) where {T<:BlasReal} = - @invoke /(B::AbstractMatrix, F::Factorization) -(/)(B::Adjoint{Complex{T},Vector{Complex{T}}}, F::AdjointFactorization{T}) where {T<:BlasReal} = - (F' \ B')' - -function (\)(F::Factorization, B::AbstractVecOrMat) - require_one_based_indexing(B) - TFB = typeof(oneunit(eltype(F)) \ oneunit(eltype(B))) - ldiv!(F, copy_similar(B, TFB)) -end -(\)(F::TransposeFactorization, B::AbstractVecOrMat) = conj!(adjoint(F.parent) \ conj.(B)) - -function (/)(B::AbstractMatrix, F::Factorization) - require_one_based_indexing(B) - TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F))) - rdiv!(copy_similar(B, TFB), F) -end -(/)(A::AbstractMatrix, F::AdjointFactorization) = adjoint(adjoint(F) \ adjoint(A)) -(/)(A::AbstractMatrix, F::TransposeFactorization) = transpose(transpose(F) \ transpose(A)) - function ldiv!(Y::AbstractVector, A::Factorization, B::AbstractVector) require_one_based_indexing(Y, B) m, n = size(A) @@ -202,3 +180,25 @@ function ldiv!(Y::AbstractMatrix, A::Factorization, B::AbstractMatrix) return ldiv!(A, Y) end end + +function (/)(B::AbstractMatrix, F::Factorization) + require_one_based_indexing(B) + TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F))) + rdiv!(copy_similar(B, TFB), F) +end +# reinterpretation trick for complex lhs and real factorization +function (/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}}}, F::Factorization{T}) where {T<:BlasReal} + require_one_based_indexing(B) + x = rdiv!(copy(reinterpret(T, B)), F) + return copy(reinterpret(Complex{T}, x)) +end +# don't do the reinterpretation for [Adjoint/Transpose]Factorization +(/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}}}, F::TransposeFactorization{T}) where {T<:BlasReal} = + @invoke /(B::AbstractMatrix, F::Factorization) +(/)(B::Matrix{Complex{T}}, F::AdjointFactorization{T}) where {T<:BlasReal} = + @invoke /(B::AbstractMatrix, F::Factorization) +(/)(B::Adjoint{Complex{T},Vector{Complex{T}}}, F::AdjointFactorization{T}) where {T<:BlasReal} = + (F' \ B')' + +rdiv!(B::AbstractMatrix, A::TransposeFactorization) = transpose(ldiv!(A.parent, transpose(B))) +rdiv!(B::AbstractMatrix, A::AdjointFactorization) = adjoint(ldiv!(A.parent, adjoint(B))) diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index 8b45a116c1a48..6f5efd1c42d16 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -600,7 +600,6 @@ function rdiv!(B::AbstractVecOrMat{<:Complex}, F::Hessenberg{<:Complex,<:Any,<:A end ldiv!(F::AdjointFactorization{<:Any,<:Hessenberg}, B::AbstractVecOrMat) = rdiv!(B', F')' -rdiv!(B::AbstractMatrix, F::AdjointFactorization{<:Any,<:Hessenberg}) = ldiv!(F', B')' det(F::Hessenberg) = det(F.H; shift=F.μ) logabsdet(F::Hessenberg) = logabsdet(F.H; shift=F.μ) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index f45386fa7649c..b01ea13ca87d3 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -709,8 +709,6 @@ function ldiv!(adjA::AdjointFactorization{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::Ab end rdiv!(B::AbstractMatrix, A::LU) = transpose(ldiv!(transpose(A), transpose(B))) -rdiv!(B::AbstractMatrix, A::TransposeFactorization{<:Any,<:LU}) = transpose(ldiv!(A.parent, transpose(B))) -rdiv!(B::AbstractMatrix, A::AdjointFactorization{<:Any,<:LU}) = adjoint(ldiv!(A.parent, adjoint(B))) # Conversions AbstractMatrix(F::LU) = (F.L * F.U)[invperm(F.p),:] From 1d3621263a88e3b5bc7105cb7095ce432dfdb0e5 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 26 Jul 2023 07:44:46 +0200 Subject: [PATCH 3/4] Apply suggestions from code review Co-authored-by: Martin Holters --- stdlib/LinearAlgebra/src/factorization.jl | 6 ++++-- stdlib/LinearAlgebra/test/hessenberg.jl | 1 + 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index bd93fccadb0ab..5b5db9ff58250 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -187,18 +187,20 @@ function (/)(B::AbstractMatrix, F::Factorization) rdiv!(copy_similar(B, TFB), F) end # reinterpretation trick for complex lhs and real factorization -function (/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}}}, F::Factorization{T}) where {T<:BlasReal} +function (/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}},Transpose{Complex{T},Vector{Complex{T}}}}, F::Factorization{T}) where {T<:BlasReal} require_one_based_indexing(B) x = rdiv!(copy(reinterpret(T, B)), F) return copy(reinterpret(Complex{T}, x)) end # don't do the reinterpretation for [Adjoint/Transpose]Factorization -(/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}}}, F::TransposeFactorization{T}) where {T<:BlasReal} = +(/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}},Transpose{Complex{T},Vector{Complex{T}}}}, F::TransposeFactorization{T}) where {T<:BlasReal} = @invoke /(B::AbstractMatrix, F::Factorization) (/)(B::Matrix{Complex{T}}, F::AdjointFactorization{T}) where {T<:BlasReal} = @invoke /(B::AbstractMatrix, F::Factorization) (/)(B::Adjoint{Complex{T},Vector{Complex{T}}}, F::AdjointFactorization{T}) where {T<:BlasReal} = (F' \ B')' +(/)(B::Transpose{Complex{T},Vector{Complex{T}}}, F::TransposeFactorization{T}) where {T<:BlasReal} = + transpose(transpose(F) \ transpose(B)) rdiv!(B::AbstractMatrix, A::TransposeFactorization) = transpose(ldiv!(A.parent, transpose(B))) rdiv!(B::AbstractMatrix, A::AdjointFactorization) = adjoint(ldiv!(A.parent, adjoint(B))) diff --git a/stdlib/LinearAlgebra/test/hessenberg.jl b/stdlib/LinearAlgebra/test/hessenberg.jl index 62d10487b9b45..105b9f8970ec8 100644 --- a/stdlib/LinearAlgebra/test/hessenberg.jl +++ b/stdlib/LinearAlgebra/test/hessenberg.jl @@ -179,6 +179,7 @@ let n = 10 @test (H - I) \ B ≈ (A - I) \ B @test (H - (3+4im)I) \ B ≈ (A - (3+4im)I) \ B @test b' / H ≈ b' / A ≈ complex(b') / H + @test transpose(b) / H ≈ transpose(b) / A ≈ transpose(complex(b)) / H @test B' / H ≈ B' / A ≈ complex(B') / H @test b' / H' ≈ complex(b)' / H' @test B' / (H - I) ≈ B' / (A - I) From 6c4b22f3d18feccab62de7129bd663d2c412607f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 26 Jul 2023 13:57:19 +0200 Subject: [PATCH 4/4] use alias (middle-ground?) --- stdlib/LinearAlgebra/src/factorization.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index 5b5db9ff58250..6f5a631cf9164 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -187,13 +187,13 @@ function (/)(B::AbstractMatrix, F::Factorization) rdiv!(copy_similar(B, TFB), F) end # reinterpretation trick for complex lhs and real factorization -function (/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}},Transpose{Complex{T},Vector{Complex{T}}}}, F::Factorization{T}) where {T<:BlasReal} +function (/)(B::Union{Matrix{Complex{T}},AdjOrTrans{Complex{T},Vector{Complex{T}}}}, F::Factorization{T}) where {T<:BlasReal} require_one_based_indexing(B) x = rdiv!(copy(reinterpret(T, B)), F) return copy(reinterpret(Complex{T}, x)) end # don't do the reinterpretation for [Adjoint/Transpose]Factorization -(/)(B::Union{Matrix{Complex{T}},Adjoint{Complex{T},Vector{Complex{T}}},Transpose{Complex{T},Vector{Complex{T}}}}, F::TransposeFactorization{T}) where {T<:BlasReal} = +(/)(B::Union{Matrix{Complex{T}},AdjOrTrans{Complex{T},Vector{Complex{T}}}}, F::TransposeFactorization{T}) where {T<:BlasReal} = @invoke /(B::AbstractMatrix, F::Factorization) (/)(B::Matrix{Complex{T}}, F::AdjointFactorization{T}) where {T<:BlasReal} = @invoke /(B::AbstractMatrix, F::Factorization)