Skip to content

Commit

Permalink
Define / method for adjoint LU rhs. (#33209)
Browse files Browse the repository at this point in the history
Fixes #33177

(cherry picked from commit 4c8f8fa)
  • Loading branch information
andreasnoack authored and KristofferC committed Sep 18, 2019
1 parent b1eefc9 commit 9000bd8
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 4 deletions.
24 changes: 20 additions & 4 deletions stdlib/LinearAlgebra/src/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -416,13 +416,29 @@ function ldiv!(adjA::Adjoint{<:Any,<:LU{<:Any,<:StridedMatrix}}, B::StridedVecOr
_apply_inverse_ipiv_rows!(A, B)
end

\(A::Adjoint{<:Any,<:LU}, B::Adjoint{<:Any,<:StridedVecOrMat}) = A \ copy(B)
\(A::Transpose{<:Any,<:LU}, B::Transpose{<:Any,<:StridedVecOrMat}) = A \ copy(B)
\(A::Adjoint{T,<:LU{T,<:StridedMatrix}}, B::Adjoint{T,<:StridedVecOrMat{T}}) where {T<:BlasComplex} =
(\)(A::Adjoint{<:Any,<:LU}, B::Adjoint{<:Any,<:StridedVecOrMat}) = A \ copy(B)
(\)(A::Transpose{<:Any,<:LU}, B::Transpose{<:Any,<:StridedVecOrMat}) = A \ copy(B)
(\)(A::Adjoint{T,<:LU{T,<:StridedMatrix}}, B::Adjoint{T,<:StridedVecOrMat{T}}) where {T<:BlasComplex} =
LAPACK.getrs!('C', A.parent.factors, A.parent.ipiv, copy(B))
\(A::Transpose{T,<:LU{T,<:StridedMatrix}}, B::Transpose{T,<:StridedVecOrMat{T}}) where {T<:BlasFloat} =
(\)(A::Transpose{T,<:LU{T,<:StridedMatrix}}, B::Transpose{T,<:StridedVecOrMat{T}}) where {T<:BlasFloat} =
LAPACK.getrs!('T', A.parent.factors, A.parent.ipiv, copy(B))

function (/)(A::AbstractMatrix, F::Adjoint{<:Any,<:LU})
T = promote_type(eltype(A), eltype(F))
return adjoint(ldiv!(F.parent, copy_oftype(adjoint(A), T)))
end
# To avoid ambiguities with definitions in adjtrans.jl and factorizations.jl
(/)(adjA::Adjoint{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:LU}) = adjoint(F.parent \ adjA.parent)
(/)(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))))
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))))
end

function det(F::LU{T}) where T
n = checksquare(F)
issuccess(F) || return zero(T)
Expand Down
16 changes: 16 additions & 0 deletions stdlib/LinearAlgebra/test/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -316,4 +316,20 @@ end
0 0 0 0 1 1 0 1]) 6
end

@testset "Issue #33177. No ldiv!(LU, Adjoint)" begin
A = [1 0; 1 1]
B = [1 2; 2 8]
F = lu(B)
@test (A / F') * B == A
@test (A' / F') * B == A'

a = complex.(randn(2), randn(2))
@test (a' / F') * B a'
@test (transpose(a) / F') * B transpose(a)

A = complex.(randn(2, 2), randn(2, 2))
@test (A' / F') * B A'
@test (transpose(A) / F') * B transpose(A)
end

end # module TestLU

0 comments on commit 9000bd8

Please sign in to comment.