diff --git a/src/diagonal.jl b/src/diagonal.jl index 9ce9783c8..498d10d2f 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -342,7 +342,12 @@ function (*)(D::Diagonal, V::AbstractVector) return D.diag .* V end +_wrapperop_reverses_products(::Adjoint, ::Diagonal) = true +_wrapperop_reverses_products(::Transpose{<:Union{Real,Complex}}, ::Diagonal{<:Union{Real,Complex}}) = true +_wrapperop_reverses_products(A, D) = false + function mul(A::AdjOrTransAbsMat, D::Diagonal) + _wrapperop_reverses_products(A, D) || return @invoke mul(A::AbstractMatrix, D::AbstractMatrix) adj = _wrapperop(A) copy(adj(adj(D) * adj(A))) end @@ -350,6 +355,7 @@ function mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Numbe @invoke mul(A::AbstractMatrix, D::AbstractMatrix) end function mul(D::Diagonal, A::AdjOrTransAbsMat) + _wrapperop_reverses_products(A, D) || return @invoke mul(D::AbstractMatrix, A::AbstractMatrix) adj = _wrapperop(A) copy(adj(adj(A) * adj(D))) end @@ -370,6 +376,7 @@ end # A' = A' * D => A = D' * A # This uses the fact that D' is a Diagonal function rmul!(A::AdjOrTransAbsMat, D::Diagonal) + _wrapperop_reverses_products(A, D) || return @invoke rmul!(A::AbstractMatrix, D::Diagonal) f = _wrapperop(A) lmul!(f(D), f(A)) A @@ -420,6 +427,7 @@ end # A' = D * A' => A = A * D' # This uses the fact that D' is a Diagonal function lmul!(D::Diagonal, A::AdjOrTransAbsMat) + _wrapperop_reverses_products(A, D) || return @invoke lmul!(D::Diagonal, A::AbstractVecOrMat) f = _wrapperop(A) rmul!(f(A), f(D)) A diff --git a/src/triangular.jl b/src/triangular.jl index 93fa0fcdc..a97525f3d 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -1841,7 +1841,7 @@ function generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, ::Function, xA: for i in axes(B,1)[2:end] C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,1])) * C1 end - for j in axes(A,2) + for j in axes(A,2)[2:end] Cj = C[j,k] for i in j+1:lastindex(A,1) C[i,k] -= _ustrip(conj(A[i,j])) * Cj diff --git a/test/diagonal.jl b/test/diagonal.jl index ef438037a..1793c63d5 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -16,6 +16,7 @@ using Main.LinearAlgebraTestHelpers.InfiniteArrays using Main.LinearAlgebraTestHelpers.FillArrays using Main.LinearAlgebraTestHelpers.SizedArrays using Main.LinearAlgebraTestHelpers.ImmutableArrays +using Main.LinearAlgebraTestHelpers.Quaternions const n=12 # Size of matrix problem to test Random.seed!(1) @@ -1578,4 +1579,17 @@ end @test D == D2 end +@testset "products of Diagonal with transposed non-commutative matrices (issue #1750)" begin + i, j, k = Quaternion(0.0,1,0,0), Quaternion(0.0,0,1,0), Quaternion(0.0,0,0,1) + M = [i j; zero(i) k] + A = view(M, [1, 2], :) + D = Diagonal([j, k]) + for op in (transpose, adjoint) + @test op(A) * D == Matrix(op(A)) * D + @test D * op(A) == D * Matrix(op(A)) + @test rmul!(op(copy(M)), D) == Matrix(op(M)) * D + @test lmul!(D, op(copy(M))) == D * Matrix(op(M)) + end +end + end # module TestDiagonal diff --git a/test/triangular.jl b/test/triangular.jl index e24f62ab6..65abf42d8 100644 --- a/test/triangular.jl +++ b/test/triangular.jl @@ -1118,4 +1118,17 @@ end end end +@testset "issue #1749: ldiv! with lazy adjoint/transpose parent" begin + M = ComplexF64[1 2+3im 1im; 4-1im 5 2-im; 3+im 2im 7] + b = ComplexF64[1, 2, 3] + B = ComplexF64[1 4; 2 5; 3 6] + for P in (transpose(M'), adjoint(transpose(M)), transpose(M), M') + for W in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular) + A = W(P) + @test A \ b ≈ Matrix(A) \ b ≈ lu(Matrix(A)) \ b + @test A \ B ≈ Matrix(A) \ B ≈ lu(Matrix(A)) \ B + end + end +end + end # module TestTriangular