diff --git a/base/essentials.jl b/base/essentials.jl index fb4608fdaa2a3..60a70f5459aee 100644 --- a/base/essentials.jl +++ b/base/essentials.jl @@ -371,7 +371,7 @@ julia> sizeof(ComplexF64) If `T` does not have a specific size, an error is thrown. ```jldoctest -julia> sizeof(Base.LinAlg.LU) +julia> sizeof(AbstractArray) ERROR: argument is an abstract type; size is indeterminate Stacktrace: [...] diff --git a/base/hashing2.jl b/base/hashing2.jl index 83ef9ea09683f..55d67ca9caef9 100644 --- a/base/hashing2.jl +++ b/base/hashing2.jl @@ -178,4 +178,4 @@ function hash(s::Union{String,SubString{String}}, h::UInt) # note: use pointer(s) here (see #6058). ccall(memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), pointer(s), sizeof(s), h % UInt32) + h end -hash(s::AbstractString, h::UInt) = hash(String(s), h) \ No newline at end of file +hash(s::AbstractString, h::UInt) = hash(String(s), h) diff --git a/base/reflection.jl b/base/reflection.jl index 1d8a725610d5f..b25f954a532bf 100644 --- a/base/reflection.jl +++ b/base/reflection.jl @@ -9,8 +9,8 @@ Get the name of a `Module` as a `Symbol`. # Examples ```jldoctest -julia> module_name(Base.LinAlg) -:LinAlg +julia> module_name(Base) +:Base ``` """ module_name(m::Module) = ccall(:jl_module_name, Ref{Symbol}, (Any,), m) @@ -25,8 +25,8 @@ Get a module's enclosing `Module`. `Main` is its own parent. julia> parentmodule(Main) Main -julia> parentmodule(Base.LinAlg.BLAS) -Base.LinAlg +julia> parentmodule(Base.Sys) +Base ``` """ parentmodule(m::Module) = ccall(:jl_module_parent, Ref{Module}, (Any,), m) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index c82d6f67d8b9a..4d18d7d5d00fd 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -240,11 +240,9 @@ function char_uplo(uplo::Symbol) end """ - ldiv!([Y,] A, B) -> Y + ldiv!(Y, A, B) -> Y Compute `A \\ B` in-place and store the result in `Y`, returning the result. -If only two arguments are passed, then `ldiv!(A, B)` overwrites `B` with -the result. The argument `A` should *not* be a matrix. Rather, instead of matrices it should be a factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)). @@ -256,11 +254,24 @@ control over the factorization of `A`. ldiv!(Y, A, B) """ - rdiv!([Y,] A, B) -> Y + ldiv!(A, B) -Compute `A / B` in-place and store the result in `Y`, returning the result. -If only two arguments are passed, then `rdiv!(A, B)` overwrites `A` with -the result. +Compute `A \\ B` in-place and overwriting `B` to store the result. + +The argument `A` should *not* be a matrix. Rather, instead of matrices it should be a +factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)). +The reason for this is that factorization itself is both expensive and typically allocates memory +(although it can also be done in-place via, e.g., [`lufact!`](@ref)), +and performance-critical situations requiring `ldiv!` usually also require fine-grained +control over the factorization of `A`. +""" +ldiv!(A, B) + + +""" + rdiv!(A, B) + +Compute `A / B` in-place and overwriting `A` to store the result. The argument `B` should *not* be a matrix. Rather, instead of matrices it should be a factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)). @@ -269,7 +280,7 @@ The reason for this is that factorization itself is both expensive and typically and performance-critical situations requiring `rdiv!` usually also require fine-grained control over the factorization of `B`. """ -rdiv!(Y, A, B) +rdiv!(A, B) copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A) copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A) diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index 6e63a0d353016..3859eb5bdfa86 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -727,11 +727,11 @@ export A_ldiv_B!, @deprecate A_mul_B!(C::StridedVecOrMat, S::SymTridiagonal, B::StridedVecOrMat) mul!(C, S, B) # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/diagonal.jl, to deprecate -@deprecate A_mul_B!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) mul!(A, D) -@deprecate A_mul_B!(A::UnitLowerTriangular, D::Diagonal) mul!(A, D) -@deprecate A_mul_B!(A::UnitUpperTriangular, D::Diagonal) mul!(A, D) -@deprecate A_mul_B!(D::Diagonal, B::UnitLowerTriangular) mul!(D, B) -@deprecate A_mul_B!(D::Diagonal, B::UnitUpperTriangular) mul!(D, B) +@deprecate A_mul_B!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) mul1!(A, D) +@deprecate A_mul_B!(A::UnitLowerTriangular, D::Diagonal) mul1!(A, D) +@deprecate A_mul_B!(A::UnitUpperTriangular, D::Diagonal) mul1!(A, D) +@deprecate A_mul_B!(D::Diagonal, B::UnitLowerTriangular) mul2!(D, B) +@deprecate A_mul_B!(D::Diagonal, B::UnitUpperTriangular) mul2!(D, B) @deprecate Ac_mul_B(D::Diagonal, B::Diagonal) (*)(adjoint(D), B) @deprecate Ac_mul_B(A::AbstractTriangular, D::Diagonal) (*)(adjoint(A), D) @deprecate Ac_mul_B(A::AbstractMatrix, D::Diagonal) (*)(adjoint(A), D) @@ -747,16 +747,19 @@ export A_ldiv_B!, @deprecate A_mul_Bt(D::Diagonal, A::AbstractMatrix) (*)(D, transpose(A)) @deprecate Ac_mul_Bc(D::Diagonal, B::Diagonal) (*)(adjoint(D), adjoint(B)) @deprecate At_mul_Bt(D::Diagonal, B::Diagonal) (*)(transpose(D), transpose(B)) -@deprecate A_mul_B!(A::Diagonal,B::Diagonal) mul!(A, B) -@deprecate At_mul_B!(A::Diagonal,B::Diagonal) mul!(transpose(A), B) -@deprecate Ac_mul_B!(A::Diagonal,B::Diagonal) mul!(adjoint(A), B) -@deprecate A_mul_B!(A::QRPackedQ, D::Diagonal) mul!(A, D) -@deprecate A_mul_B!(A::Diagonal,B::AbstractMatrix) mul!(A, B) -@deprecate At_mul_B!(A::Diagonal,B::AbstractMatrix) mul!(transpose(A), B) -@deprecate Ac_mul_B!(A::Diagonal,B::AbstractMatrix) mul!(adjoint(A), B) -@deprecate A_mul_B!(A::AbstractMatrix,B::Diagonal) mul!(A, B) -@deprecate A_mul_Bt!(A::AbstractMatrix,B::Diagonal) mul!(A, transpose(B)) -@deprecate A_mul_Bc!(A::AbstractMatrix,B::Diagonal) mul!(A, adjoint(B)) +function A_mul_B!(A::Diagonal,B::Diagonal) + depwarn("`A_mul_B!(A::Diagonal,B::Diagonal)` should be replaced with `mul1!(A, B)` or `mul2!(A, B)`.", :A_mul_B!) + throw(MethodError(A_mul_B!, (A, B))) +end +@deprecate At_mul_B!(A::Diagonal,B::Diagonal) mul2!(transpose(A), B) +@deprecate Ac_mul_B!(A::Diagonal,B::Diagonal) mul2!(adjoint(A), B) +@deprecate A_mul_B!(A::QRPackedQ, D::Diagonal) mul1!(A, D) +@deprecate A_mul_B!(A::Diagonal,B::AbstractMatrix) mul2!(A, B) +@deprecate At_mul_B!(A::Diagonal,B::AbstractMatrix) mul2!(transpose(A), B) +@deprecate Ac_mul_B!(A::Diagonal,B::AbstractMatrix) mul2!(adjoint(A), B) +@deprecate A_mul_B!(A::AbstractMatrix,B::Diagonal) mul1!(A, B) +@deprecate A_mul_Bt!(A::AbstractMatrix,B::Diagonal) mul1!(A, transpose(B)) +@deprecate A_mul_Bc!(A::AbstractMatrix,B::Diagonal) mul1!(A, adjoint(B)) @deprecate A_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector) mul!(out, A, in) @deprecate Ac_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector) mul!(out, adjoint(A), in) @deprecate At_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector) mul!(out, transpose(A), in) @@ -780,7 +783,7 @@ export A_ldiv_B!, @deprecate A_ldiv_B!(D::Diagonal, B::StridedVecOrMat) ldiv!(D, B) # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/special.jl, to deprecate -@deprecate A_mul_Bc!(A::AbstractTriangular, B::Union{QRCompactWYQ,QRPackedQ}) mul!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::AbstractTriangular, B::Union{QRCompactWYQ,QRPackedQ}) mul1!(A, adjoint(B)) @deprecate A_mul_Bc(A::AbstractTriangular, B::Union{QRCompactWYQ,QRPackedQ}) (*)(A, adjoint(B)) # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/bunchkaufman.jl, to deprecate @@ -805,10 +808,10 @@ export A_ldiv_B!, @deprecate At_ldiv_B(F::Factorization, B) (\)(transpose(F), B) # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/hessenberg.jl, to deprecate -@deprecate A_mul_B!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} mul!(Q, X) -@deprecate A_mul_B!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} mul!(X, Q) -@deprecate Ac_mul_B!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} mul!(adjoint(Q), X) -@deprecate A_mul_Bc!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} mul!(X, adjoint(Q)) +@deprecate A_mul_B!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} mul2!(Q, X) +@deprecate A_mul_B!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} mul1!(X, Q) +@deprecate Ac_mul_B!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} mul2!(adjoint(Q), X) +@deprecate A_mul_Bc!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} mul1!(X, adjoint(Q)) @deprecate Ac_mul_B(Q::HessenbergQ{T}, X::StridedVecOrMat{S}) where {T,S} (*)(adjoint(Q), X) @deprecate A_mul_Bc(X::StridedVecOrMat{S}, Q::HessenbergQ{T}) where {T,S} (*)(X, adjoint(Q)) @@ -858,43 +861,43 @@ export A_ldiv_B!, @deprecate Ac_ldiv_B!(A::LU{T,Tridiagonal{T,V}}, B::AbstractVecOrMat) where {T,V} ldiv!(adjoint(A), B) # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/lq.jl, to deprecate -@deprecate A_mul_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} mul!(A, B) -@deprecate A_mul_B!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat} mul!(A, B) -@deprecate A_mul_B!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat} mul!(A, B) -@deprecate A_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} mul!(A, B) -@deprecate Ac_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasReal} mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasComplex} mul!(adjoint(A), B) +@deprecate A_mul_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} mul2!(A, B) +@deprecate A_mul_B!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat} A*B +@deprecate A_mul_B!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat} A*B +@deprecate A_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} mul2!(A, B) +@deprecate Ac_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasReal} mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasComplex} mul2!(adjoint(A), B) @deprecate Ac_mul_B(A::LQPackedQ, B::StridedVecOrMat) (*)(adjoint(A), B) @deprecate A_mul_Bc(A::LQPackedQ, B::StridedVecOrMat) (*)(A, adjoint(B)) @deprecate Ac_mul_Bc(A::LQPackedQ, B::StridedVecOrMat) (*)(adjoint(A), adjoint(B)) -@deprecate A_mul_B!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} mul!(A, B) -@deprecate A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasReal} mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasComplex} mul!(A, adjoint(B)) +@deprecate A_mul_B!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} mul1!(A, B) +@deprecate A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasReal} mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasComplex} mul1!(A, adjoint(B)) @deprecate A_mul_Bc(A::StridedVecOrMat, Q::LQPackedQ) (*)(A, adjoint(Q)) @deprecate Ac_mul_Bc(A::StridedMatrix, Q::LQPackedQ) (*)(adjoint(A), adjoint(Q)) @deprecate Ac_mul_B(A::StridedMatrix, Q::LQPackedQ) (*)(adjoint(A), Q) @deprecate A_ldiv_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T} ldiv!(A, B) # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/qr.jl, to deprecate -@deprecate A_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} mul!(A, B) -@deprecate A_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} mul!(A, B) -@deprecate A_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) mul!(A, B) -@deprecate Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) mul!(adjoint(A), B) +@deprecate A_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} mul2!(A, B) +@deprecate A_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} mul2!(A, B) +@deprecate A_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) mul2!(A, B) +@deprecate Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) mul2!(adjoint(A), B) @deprecate Ac_mul_B(Q::AbstractQ, B::StridedVecOrMat) (*)(adjoint(Q), B) @deprecate A_mul_Bc(Q::AbstractQ, B::StridedVecOrMat) (*)(Q, adjoint(B)) @deprecate Ac_mul_Bc(Q::AbstractQ, B::StridedVecOrMat) (*)(adjoint(Q), adjoint(B)) -@deprecate A_mul_B!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} mul!(A, B) -@deprecate A_mul_B!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} mul!(A, B) -@deprecate A_mul_B!(A::StridedMatrix,Q::QRPackedQ) mul!(A, Q) -@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasReal} mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasComplex} mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasReal} mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasComplex} mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::StridedMatrix,Q::QRPackedQ) mul!(A, adjoint(Q)) +@deprecate A_mul_B!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} mul1!(A, B) +@deprecate A_mul_B!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} mul1!(A, B) +@deprecate A_mul_B!(A::StridedMatrix,Q::QRPackedQ) mul1!(A, Q) +@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasReal} mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasComplex} mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasReal} mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasComplex} mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::StridedMatrix,Q::QRPackedQ) mul1!(A, adjoint(Q)) @deprecate A_mul_Bc(A::StridedMatrix, B::AbstractQ) (*)(A, adjoint(B)) @deprecate A_mul_Bc(rowvec::RowVector, B::AbstractQ) (*)(rowvec, adjoint(B)) @deprecate Ac_mul_B(A::StridedVecOrMat, Q::AbstractQ) (*)(adjoint(A), Q) @@ -987,20 +990,20 @@ export A_ldiv_B!, @deprecate At_mul_Bt(A::AbstractTriangular, B::AbstractTriangular) (*)(transpose(A), transpose(B)) @deprecate At_mul_Bt(A::AbstractTriangular, B::AbstractMatrix) (*)(transpose(A), transpose(B)) @deprecate At_mul_Bt(A::AbstractMatrix, B::AbstractTriangular) (*)(transpose(A), transpose(B)) -@deprecate A_mul_Bc!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) mul!(A, adjoint(B)) -@deprecate A_mul_Bt!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) mul!(A, transpose(B)) -@deprecate A_mul_Bt!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) mul!(A, transpose(B)) +@deprecate A_mul_Bc!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) mul1!(A, adjoint(B)) +@deprecate A_mul_Bt!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) mul1!(A, transpose(B)) +@deprecate A_mul_Bt!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) mul1!(A, transpose(B)) @deprecate A_rdiv_Bc!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) rdiv!(A, adjoint(B)) @deprecate A_rdiv_Bc!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) rdiv!(A, adjoint(B)) @deprecate A_rdiv_Bt!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) rdiv!(A, transpose(B)) @deprecate A_rdiv_Bt!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) rdiv!(A, transpose(B)) @deprecate A_rdiv_B!(A::UpperTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) rdiv!(A, B) @deprecate A_rdiv_B!(A::LowerTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) rdiv!(A, B) -@deprecate Ac_mul_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular) mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular) mul!(adjoint(A), B) -@deprecate At_mul_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular) mul!(transpose(A), B) -@deprecate At_mul_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular) mul!(transpose(A), B) +@deprecate Ac_mul_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular) mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular) mul2!(adjoint(A), B) +@deprecate At_mul_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular) mul2!(transpose(A), B) +@deprecate At_mul_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular) mul2!(transpose(A), B) @deprecate Ac_ldiv_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular) ldiv!(adjoint(A), B) @deprecate Ac_ldiv_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular) ldiv!(adjoint(A), B) @deprecate At_ldiv_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular) ldiv!(transpose(A), B) @@ -1025,30 +1028,30 @@ export A_ldiv_B!, @deprecate At_ldiv_B!(A::UpperTriangular, b::AbstractVector, x::AbstractVector = b) ldiv!(transpose(A), b, x) @deprecate At_ldiv_B!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector = b) ldiv!(transpose(A), b, x) @deprecate At_ldiv_B!(A::LowerTriangular, b::AbstractVector, x::AbstractVector = b) ldiv!(transpose(A), b, x) -@deprecate A_mul_Bt!(A::StridedMatrix, B::UnitLowerTriangular) mul!(A, transpose(B)) -@deprecate A_mul_Bt!(A::StridedMatrix, B::LowerTriangular) mul!(A, transpose(B)) -@deprecate A_mul_Bt!(A::StridedMatrix, B::UnitUpperTriangular) mul!(A, transpose(B)) -@deprecate A_mul_Bt!(A::StridedMatrix, B::UpperTriangular) mul!(A, transpose(B)) -@deprecate A_mul_Bc!(A::StridedMatrix, B::UnitLowerTriangular) mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::StridedMatrix, B::LowerTriangular) mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::StridedMatrix, B::UnitUpperTriangular) mul!(A, adjoint(B)) -@deprecate A_mul_Bc!(A::StridedMatrix, B::UpperTriangular) mul!(A, adjoint(B)) -@deprecate A_mul_B!(A::StridedMatrix, B::UnitLowerTriangular) mul!(A, B) -@deprecate A_mul_B!(A::StridedMatrix, B::LowerTriangular) mul!(A, B) -@deprecate A_mul_B!(A::StridedMatrix, B::UnitUpperTriangular) mul!(A, B) -@deprecate A_mul_B!(A::StridedMatrix, B::UpperTriangular) mul!(A, B) -@deprecate At_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) mul!(transpose(A), B) -@deprecate At_mul_B!(A::LowerTriangular, B::StridedVecOrMat) mul!(transpose(A), B) -@deprecate At_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) mul!(transpose(A), B) -@deprecate At_mul_B!(A::UpperTriangular, B::StridedVecOrMat) mul!(transpose(A), B) -@deprecate Ac_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::LowerTriangular, B::StridedVecOrMat) mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) mul!(adjoint(A), B) -@deprecate Ac_mul_B!(A::UpperTriangular, B::StridedVecOrMat) mul!(adjoint(A), B) -@deprecate A_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) mul!(A, B) -@deprecate A_mul_B!(A::LowerTriangular, B::StridedVecOrMat) mul!(A, B) -@deprecate A_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) mul!(A, B) -@deprecate A_mul_B!(A::UpperTriangular, B::StridedVecOrMat) mul!(A, B) +@deprecate A_mul_Bt!(A::StridedMatrix, B::UnitLowerTriangular) mul1!(A, transpose(B)) +@deprecate A_mul_Bt!(A::StridedMatrix, B::LowerTriangular) mul1!(A, transpose(B)) +@deprecate A_mul_Bt!(A::StridedMatrix, B::UnitUpperTriangular) mul1!(A, transpose(B)) +@deprecate A_mul_Bt!(A::StridedMatrix, B::UpperTriangular) mul1!(A, transpose(B)) +@deprecate A_mul_Bc!(A::StridedMatrix, B::UnitLowerTriangular) mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::StridedMatrix, B::LowerTriangular) mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::StridedMatrix, B::UnitUpperTriangular) mul1!(A, adjoint(B)) +@deprecate A_mul_Bc!(A::StridedMatrix, B::UpperTriangular) mul1!(A, adjoint(B)) +@deprecate A_mul_B!(A::StridedMatrix, B::UnitLowerTriangular) mul1!(A, B) +@deprecate A_mul_B!(A::StridedMatrix, B::LowerTriangular) mul1!(A, B) +@deprecate A_mul_B!(A::StridedMatrix, B::UnitUpperTriangular) mul1!(A, B) +@deprecate A_mul_B!(A::StridedMatrix, B::UpperTriangular) mul1!(A, B) +@deprecate At_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) mul2!(transpose(A), B) +@deprecate At_mul_B!(A::LowerTriangular, B::StridedVecOrMat) mul2!(transpose(A), B) +@deprecate At_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) mul2!(transpose(A), B) +@deprecate At_mul_B!(A::UpperTriangular, B::StridedVecOrMat) mul2!(transpose(A), B) +@deprecate Ac_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::LowerTriangular, B::StridedVecOrMat) mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) mul2!(adjoint(A), B) +@deprecate Ac_mul_B!(A::UpperTriangular, B::StridedVecOrMat) mul2!(adjoint(A), B) +@deprecate A_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) mul2!(A, B) +@deprecate A_mul_B!(A::LowerTriangular, B::StridedVecOrMat) mul2!(A, B) +@deprecate A_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) mul2!(A, B) +@deprecate A_mul_B!(A::UpperTriangular, B::StridedVecOrMat) mul2!(A, B) @deprecate A_mul_B!(C::AbstractVector , A::AbstractTriangular, B::AbstractVector) mul!(C, A, B) @deprecate A_mul_B!(C::AbstractMatrix , A::AbstractTriangular, B::AbstractVecOrMat) mul!(C, A, B) @deprecate A_mul_B!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) mul!(C, A, B) @@ -1058,7 +1061,7 @@ export A_ldiv_B!, @deprecate At_mul_B!(C::AbstractVector , A::AbstractTriangular, B::AbstractVector) mul!(C, transpose(A), B) @deprecate At_mul_B!(C::AbstractMatrix , A::AbstractTriangular, B::AbstractVecOrMat) mul!(C, transpose(A), B) @deprecate At_mul_B!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) mul!(C, transpose(A), B) -@deprecate A_mul_B!(A::Tridiagonal, B::AbstractTriangular) mul!(A, B) +@deprecate A_mul_B!(A::Tridiagonal, B::AbstractTriangular) mul2!(A, B) @deprecate A_mul_B!(C::AbstractMatrix, A::AbstractTriangular, B::Tridiagonal) mul!(C, A, B) @deprecate A_mul_B!(C::AbstractMatrix, A::Tridiagonal, B::AbstractTriangular) mul!(C, A, B) @deprecate A_mul_Bt!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) mul!(C, A, transpose(B)) @@ -1110,22 +1113,22 @@ for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'), (:UnitUpperTriangular, 'U', 'U')) @eval begin # Vector multiplication - @deprecate A_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} mul!(A, b) - @deprecate At_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} mul!(transpose(A), b) - @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasReal} mul!(adjoint(A), b) - @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasComplex} mul!(adjoint(A), b) + @deprecate A_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} mul2!(A, b) + @deprecate At_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} mul2!(transpose(A), b) + @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasReal} mul2!(adjoint(A), b) + @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasComplex} mul2!(adjoint(A), b) # Matrix multiplication - @deprecate A_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} mul!(A, B) - @deprecate A_mul_B!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} mul!(A, B) + @deprecate A_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} mul2!(A, B) + @deprecate A_mul_B!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} mul2!(A, B) - @deprecate At_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} mul!(transpose(A), B) - @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasComplex} mul!(adjoint(A), B) - @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasReal} mul!(adjoint(A), B) + @deprecate At_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} mul2!(transpose(A), B) + @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasComplex} mul2!(adjoint(A), B) + @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasReal} mul2!(adjoint(A), B) - @deprecate A_mul_Bt!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} mul!(A, transpose(B)) - @deprecate A_mul_Bc!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasComplex} mul!(A, adjoint(B)) - @deprecate A_mul_Bc!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasReal} mul!(A, adjoint(B)) + @deprecate A_mul_Bt!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} mul1!(A, transpose(B)) + @deprecate A_mul_Bc!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasComplex} mul1!(A, adjoint(B)) + @deprecate A_mul_Bc!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasReal} mul1!(A, adjoint(B)) # Left division @deprecate A_ldiv_B!(A::$t{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} ldiv!(A, B) @@ -1183,12 +1186,12 @@ end @deprecate A_mul_Bt(mat::AbstractMatrix, rowvec::RowVector) (*)(mat, transpose(rowvec)) # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/givens.jl, to deprecate -@deprecate A_mul_Bc!(A::AbstractMatrix, R::Rotation) mul!(A, adjoint(R)) -@deprecate A_mul_B!(R::Rotation, A::AbstractMatrix) mul!(R, A) -@deprecate A_mul_B!(G::Givens, R::Rotation) mul!(G, R) -@deprecate A_mul_Bc!(A::AbstractMatrix, G::Givens) mul!(A, adjoint(G)) -@deprecate A_mul_B!(G::Givens, A::AbstractVecOrMat) mul!(G, A) -@deprecate A_mul_B!(G1::Givens, G2::Givens) mul!(G1, G2) +@deprecate A_mul_Bc!(A::AbstractMatrix, R::Rotation) mul1!(A, adjoint(R)) +@deprecate A_mul_B!(R::Rotation, A::AbstractMatrix) mul2!(R, A) +@deprecate A_mul_B!(G::Givens, R::Rotation) mul2!(G, R) +@deprecate A_mul_Bc!(A::AbstractMatrix, G::Givens) mul1!(A, adjoint(G)) +@deprecate A_mul_B!(G::Givens, A::AbstractVecOrMat) mul2!(G, A) +@deprecate A_mul_B!(G1::Givens, G2::Givens) G1 * G2 @deprecate A_mul_Bc(A::AbstractVecOrMat{T}, R::AbstractRotation{S}) where {T,S} (*)(A, adjoint(R)) @@ -1254,7 +1257,10 @@ norm(tv::RowVector) = norm(rvtranspose(tv)) *(A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:RowVector}) = A * rvtranspose(B.parent) *(A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:RowVector}) = A * rvadjoint(B.parent) - +@deprecate *(A::LQ,B::QR) A*Matrix(B) +@deprecate *(A::QR,B::LQ) A*Matrix(B) +@deprecate *(A::Adjoint{<:Any,<:LQ}, B::LQ) A*Matrix(B) +@deprecate *(A::LQ, B::Adjoint{<:Any,<:LQ}) A*Matrix(B) # PR #25184. Use getproperty instead of getindex for Factorizations function getindex(F::Factorization, s::Symbol) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 8e3f06e1000d1..e90820b949eed 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -151,38 +151,40 @@ end (*)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .* Db.diag) (*)(D::Diagonal, V::AbstractVector) = D.diag .* V -(*)(A::AbstractTriangular, D::Diagonal) = mul!(copy(A), D) -(*)(D::Diagonal, B::AbstractTriangular) = mul!(D, copy(B)) +(*)(A::AbstractTriangular, D::Diagonal) = mul1!(copy(A), D) +(*)(D::Diagonal, B::AbstractTriangular) = mul2!(D, copy(B)) (*)(A::AbstractMatrix, D::Diagonal) = scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A, D.diag) (*)(D::Diagonal, A::AbstractMatrix) = scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), D.diag, A) -mul!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) = typeof(A)(mul!(A.data, D)) -function mul!(A::UnitLowerTriangular, D::Diagonal) - mul!(A.data, D) + +mul1!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) = typeof(A)(mul1!(A.data, D)) +function mul1!(A::UnitLowerTriangular, D::Diagonal) + mul1!(A.data, D) for i = 1:size(A, 1) A.data[i,i] = D.diag[i] end LowerTriangular(A.data) end -function mul!(A::UnitUpperTriangular, D::Diagonal) - mul!(A.data, D) +function mul1!(A::UnitUpperTriangular, D::Diagonal) + mul1!(A.data, D) for i = 1:size(A, 1) A.data[i,i] = D.diag[i] end UpperTriangular(A.data) end -function mul!(D::Diagonal, B::UnitLowerTriangular) - mul!(D, B.data) + +function mul2!(D::Diagonal, B::UnitLowerTriangular) + mul2!(D, B.data) for i = 1:size(B, 1) B.data[i,i] = D.diag[i] end LowerTriangular(B.data) end -function mul!(D::Diagonal, B::UnitUpperTriangular) - mul!(D, B.data) +function mul2!(D::Diagonal, B::UnitUpperTriangular) + mul2!(D, B.data) for i = 1:size(B, 1) B.data[i,i] = D.diag[i] end @@ -190,40 +192,40 @@ function mul!(D::Diagonal, B::UnitUpperTriangular) end *(D::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(adjoint.(D.parent.diag) .* B.diag) -*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) = mul!(copy(A), D) +*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) = mul1!(copy(A), D) function *(adjA::Adjoint{<:Any,<:AbstractMatrix}, D::Diagonal) A = adjA.parent Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1))) adjoint!(Ac, A) - mul!(Ac, D) + mul1!(Ac, D) end *(D::Transpose{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(transpose.(D.parent.diag) .* B.diag) -*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) = mul!(copy(A), D) +*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) = mul1!(copy(A), D) function *(transA::Transpose{<:Any,<:AbstractMatrix}, D::Diagonal) A = transA.parent At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1))) transpose!(At, A) - mul!(At, D) + mul1!(At, D) end *(D::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = Diagonal(D.diag .* adjoint.(B.parent.diag)) -*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) = mul!(D, collect(B)) -*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; mul!(Array(D), adjoint(Q))) +*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) = mul2!(D, collect(B)) +*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; mul1!(Array(D), adjoint(Q))) function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix}) A = adjA.parent Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1))) adjoint!(Ac, A) - mul!(D, Ac) + mul2!(D, Ac) end *(D::Diagonal, B::Transpose{<:Any,<:Diagonal}) = Diagonal(D.diag .* transpose.(B.parent.diag)) -*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) = mul!(D, copy(B)) +*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) = mul2!(D, copy(B)) function *(D::Diagonal, transA::Transpose{<:Any,<:AbstractMatrix}) A = transA.parent At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1))) transpose!(At, A) - mul!(D, At) + mul2!(D, At) end *(D::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) = @@ -231,27 +233,15 @@ end *(D::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) = Diagonal(transpose.(D.parent.diag) .* transpose.(B.parent.diag)) -mul!(A::Diagonal, B::Diagonal) = throw(MethodError(mul!, (A, B))) -mul!(A::QRPackedQ, D::Diagonal) = throw(MethodError(mul!, (A, D))) -mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::QRPackedQ, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Diagonal) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Diagonal, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B))) -mul!(A::Transpose{<:Any,<:Diagonal}, B::Diagonal) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = throw(MethodError(mul!, (A, B))) -mul!(A::Diagonal, B::AbstractMatrix) = scale!(A.diag, B) -mul!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = adjA.parent; scale!(conj(A.diag),B)) -mul!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = transA.parent; scale!(A.diag,B)) -mul!(A::AbstractMatrix, B::Diagonal) = scale!(A,B.diag) -mul!(A::AbstractMatrix, adjB::Adjoint{<:Any,<:Diagonal}) = (B = adjB.parent; scale!(A,conj(B.diag))) -mul!(A::AbstractMatrix, transB::Transpose{<:Any,<:Diagonal}) = (B = transB.parent; scale!(A,B.diag)) +mul1!(A::Diagonal,B::Diagonal) = Diagonal(A.diag .*= B.diag) +mul2!(A::Diagonal,B::Diagonal) = Diagonal(B.diag .= A.diag .* B.diag) + +mul2!(A::Diagonal, B::AbstractMatrix) = scale!(A.diag, B) +mul2!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = adjA.parent; scale!(conj(A.diag),B)) +mul2!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = transA.parent; scale!(A.diag,B)) +mul1!(A::AbstractMatrix, B::Diagonal) = scale!(A,B.diag) +mul1!(A::AbstractMatrix, adjB::Adjoint{<:Any,<:Diagonal}) = (B = adjB.parent; scale!(A,conj(B.diag))) +mul1!(A::AbstractMatrix, transB::Transpose{<:Any,<:Diagonal}) = (B = transB.parent; scale!(A,B.diag)) # Get ambiguous method if try to unify AbstractVector/AbstractMatrix here using AbstractVecOrMat mul!(out::AbstractVector, A::Diagonal, in::AbstractVector) = out .= A.diag .* in @@ -287,6 +277,7 @@ mul!(C::AbstractMatrix, A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Rea (/)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag ./ Db.diag) + function ldiv!(D::Diagonal{T}, v::AbstractVector{T}) where {T} if length(v) != length(D.diag) throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(length(v)) rows")) @@ -316,6 +307,7 @@ function ldiv!(D::Diagonal{T}, V::AbstractMatrix{T}) where {T} V end + ldiv!(adjD::Adjoint{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} = (D = adjD.parent; ldiv!(conj(D), B)) ldiv!(transD::Transpose{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} = @@ -339,6 +331,7 @@ function rdiv!(A::AbstractMatrix{T}, D::Diagonal{T}) where {T} A end + rdiv!(A::AbstractMatrix{T}, adjD::Adjoint{<:Any,<:Diagonal{T}}) where {T} = (D = adjD.parent; rdiv!(A, conj(D))) rdiv!(A::AbstractMatrix{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} = diff --git a/stdlib/LinearAlgebra/src/givens.jl b/stdlib/LinearAlgebra/src/givens.jl index 4c84aac293974..9b508d5e1677c 100644 --- a/stdlib/LinearAlgebra/src/givens.jl +++ b/stdlib/LinearAlgebra/src/givens.jl @@ -7,14 +7,14 @@ transpose(R::AbstractRotation) = error("transpose not implemented for $(typeof(R function *(R::AbstractRotation{T}, A::AbstractVecOrMat{S}) where {T,S} TS = typeof(zero(T)*zero(S) + zero(T)*zero(S)) - mul!(convert(AbstractRotation{TS}, R), TS == S ? copy(A) : convert(AbstractArray{TS}, A)) + mul2!(convert(AbstractRotation{TS}, R), TS == S ? copy(A) : convert(AbstractArray{TS}, A)) end *(A::AbstractVector, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR) *(A::AbstractMatrix, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR) function _absvecormat_mul_adjrot(A::AbstractVecOrMat{T}, adjR::Adjoint{<:Any,<:AbstractRotation{S}}) where {T,S} R = adjR.parent TS = typeof(zero(T)*zero(S) + zero(T)*zero(S)) - mul!(TS == T ? copy(A) : convert(AbstractArray{TS}, A), adjoint(convert(AbstractRotation{TS}, R))) + mul1!(TS == T ? copy(A) : convert(AbstractArray{TS}, A), adjoint(convert(AbstractRotation{TS}, R))) end """ LinearAlgebra.Givens(i1,i2,c,s) -> G @@ -325,10 +325,7 @@ function getindex(G::Givens, i::Integer, j::Integer) end end - -mul!(G1::Givens, G2::Givens) = error("Operation not supported. Consider *") - -function mul!(G::Givens, A::AbstractVecOrMat) +function mul2!(G::Givens, A::AbstractVecOrMat) m, n = size(A, 1), size(A, 2) if G.i2 > m throw(DimensionMismatch("column indices for rotation are outside the matrix")) @@ -340,7 +337,7 @@ function mul!(G::Givens, A::AbstractVecOrMat) end return A end -function mul!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens}) +function mul1!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens}) G = adjG.parent m, n = size(A, 1), size(A, 2) if G.i2 > n @@ -353,20 +350,21 @@ function mul!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens}) end return A end -function mul!(G::Givens, R::Rotation) + +function mul2!(G::Givens, R::Rotation) push!(R.rotations, G) return R end -function mul!(R::Rotation, A::AbstractMatrix) +function mul2!(R::Rotation, A::AbstractMatrix) @inbounds for i = 1:length(R.rotations) - mul!(R.rotations[i], A) + mul2!(R.rotations[i], A) end return A end -function mul!(A::AbstractMatrix, adjR::Adjoint{<:Any,<:Rotation}) +function mul1!(A::AbstractMatrix, adjR::Adjoint{<:Any,<:Rotation}) R = adjR.parent @inbounds for i = 1:length(R.rotations) - mul!(A, adjoint(R.rotations[i])) + mul1!(A, adjoint(R.rotations[i])) end return A end @@ -388,14 +386,3 @@ end # dismabiguation methods: *(Diag/AbsTri, Adj of AbstractRotation) *(A::Diagonal, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B) *(A::AbstractTriangular, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B) -# moar disambiguation -mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B))) -mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B))) -mul!(A::Diagonal, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B))) -mul!(A::Diagonal, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B))) -mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B))) -mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B))) -mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B))) diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index 0d2c0795efd08..18a098f1fb8f7 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -73,7 +73,7 @@ function getindex(A::HessenbergQ, i::Integer, j::Integer) x[i] = 1 y = zeros(eltype(A), size(A, 2)) y[j] = 1 - return dot(x, mul!(A, y)) + return dot(x, mul2!(A, y)) end ## reconstruct the original matrix @@ -84,31 +84,30 @@ AbstractArray(F::Hessenberg) = AbstractMatrix(F) Matrix(F::Hessenberg) = Array(AbstractArray(F)) Array(F::Hessenberg) = Matrix(F) -mul!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} = +mul2!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormhr!('L', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X) -mul!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} = +mul1!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} = LAPACK.ormhr!('R', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X) -mul!(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} = +mul2!(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} = (Q = adjQ.parent; LAPACK.ormhr!('L', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X)) -mul!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T<:BlasFloat} = +mul1!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T<:BlasFloat} = (Q = adjQ.parent; LAPACK.ormhr!('R', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X)) - function (*)(Q::HessenbergQ{T}, X::StridedVecOrMat{S}) where {T,S} TT = typeof(zero(T)*zero(S) + zero(T)*zero(S)) - return mul!(Q, copy_oftype(X, TT)) + return mul2!(Q, copy_oftype(X, TT)) end function (*)(X::StridedVecOrMat{S}, Q::HessenbergQ{T}) where {T,S} TT = typeof(zero(T)*zero(S) + zero(T)*zero(S)) - return mul!(copy_oftype(X, TT), Q) + return mul1!(copy_oftype(X, TT), Q) end function *(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{S}) where {T,S} Q = adjQ.parent TT = typeof(zero(T)*zero(S) + zero(T)*zero(S)) - return mul!(adjoint(Q), copy_oftype(X, TT)) + return mul2!(adjoint(Q), copy_oftype(X, TT)) end function *(X::StridedVecOrMat{S}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T,S} Q = adjQ.parent TT = typeof(zero(T)*zero(S) + zero(T)*zero(S)) - return mul!(copy_oftype(X, TT), adjoint(Q)) + return mul1!(copy_oftype(X, TT), adjoint(Q)) end diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl index 254edae9d6036..7197c3cfb9d8b 100644 --- a/stdlib/LinearAlgebra/src/lq.jl +++ b/stdlib/LinearAlgebra/src/lq.jl @@ -57,7 +57,7 @@ function lq(A::Union{Number,AbstractMatrix}; full::Bool = false, thin::Union{Boo end F = lqfact(A) L, Q = F.L, F.Q - return L, !full ? Array(Q) : mul!(Q, Matrix{eltype(Q)}(I, size(Q.factors, 2), size(Q.factors, 2))) + return L, !full ? Array(Q) : mul2!(Q, Matrix{eltype(Q)}(I, size(Q.factors, 2), size(Q.factors, 2))) end copy(A::LQ) = LQ(copy(A.factors), copy(A.τ)) @@ -88,7 +88,7 @@ end Base.propertynames(F::LQ, private::Bool=false) = append!([:L,:Q], private ? fieldnames(typeof(F)) : Symbol[]) getindex(A::LQPackedQ, i::Integer, j::Integer) = - mul!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i] + mul2!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i] function show(io::IO, C::LQ) println(io, "$(typeof(C)) with factors L and Q:") @@ -122,47 +122,34 @@ end ## Multiplication by LQ -mul!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = - A.L * LAPACK.ormlq!('L', 'N', A.factors, A.τ, B) -mul!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat} = - A.L * LAPACK.ormlq!('L', 'N', A.factors, A.τ, Matrix(B)) -mul!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat} = - mul!(zeros(eltype(A), size(A)), Matrix(A), Matrix(B)) +mul2!(A::LQ, B::StridedVecOrMat) = + mul2!(LowerTriangular(A.L), mul2!(A.Q, B)) function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB} TAB = promote_type(TA, TB) - mul!(Factorization{TAB}(A), copy_oftype(B, TAB)) + mul2!(Factorization{TAB}(A), copy_oftype(B, TAB)) end -function *(A::LQ{TA},B::QR{TB}) where {TA,TB} - TAB = promote_type(TA, TB) - mul!(Factorization{TAB}(A), Factorization{TAB}(B)) -end -function *(A::QR{TA},B::LQ{TB}) where {TA,TB} - TAB = promote_type(TA, TB) - mul!(Factorization{TAB}(A), Factorization{TAB}(B)) -end -*(A::Adjoint{<:Any,<:LQ}, B::LQ) = copy(A) * B -*(A::LQ, B::Adjoint{<:Any,<:LQ}) = A * copy(B) ## Multiplication by Q ### QB -mul!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B) +mul2!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B) function (*)(A::LQPackedQ, B::StridedVecOrMat) TAB = promote_type(eltype(A), eltype(B)) - mul!(AbstractMatrix{TAB}(A), copy_oftype(B, TAB)) + mul2!(AbstractMatrix{TAB}(A), copy_oftype(B, TAB)) end ### QcB -mul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} = +mul2!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} = (A = adjA.parent; LAPACK.ormlq!('L','T',A.factors,A.τ,B)) -mul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = +mul2!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = (A = adjA.parent; LAPACK.ormlq!('L','C',A.factors,A.τ,B)) + function *(adjA::Adjoint{<:Any,<:LQPackedQ}, B::StridedVecOrMat) A = adjA.parent TAB = promote_type(eltype(A), eltype(B)) if size(B,1) == size(A.factors,2) - mul!(adjoint(AbstractMatrix{TAB}(A)), copy_oftype(B, TAB)) + mul2!(adjoint(AbstractMatrix{TAB}(A)), copy_oftype(B, TAB)) elseif size(B,1) == size(A.factors,1) - mul!(adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))]) + mul2!(adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))]) else throw(DimensionMismatch("first dimension of B, $(size(B,1)), must equal one of the dimensions of A, $(size(A))")) end @@ -174,14 +161,14 @@ function *(A::LQPackedQ, adjB::Adjoint{<:Any,<:StridedVecOrMat}) TAB = promote_type(eltype(A), eltype(B)) BB = similar(B, TAB, (size(B, 2), size(B, 1))) adjoint!(BB, B) - return mul!(A, BB) + return mul2!(A, BB) end function *(adjA::Adjoint{<:Any,<:LQPackedQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat}) A, B = adjA.parent, adjB.parent TAB = promote_type(eltype(A), eltype(B)) BB = similar(B, TAB, (size(B, 2), size(B, 1))) adjoint!(BB, B) - return mul!(adjoint(A), BB) + return mul2!(adjoint(A), BB) end # in-place right-application of LQPackedQs @@ -189,11 +176,11 @@ end # match the number of columns (nQ) of the LQPackedQ (Q) (necessary for in-place # operation, and the underlying LAPACK routine (ormlq) treats the implicit Q # as its (nQ-by-nQ) square form) -mul!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} = +mul1!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} = LAPACK.ormlq!('R', 'N', B.factors, B.τ, A) -mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasReal} = +mul1!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasReal} = (B = adjB.parent; LAPACK.ormlq!('R', 'T', B.factors, B.τ, A)) -mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasComplex} = +mul1!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasComplex} = (B = adjB.parent; LAPACK.ormlq!('R', 'C', B.factors, B.τ, A)) # out-of-place right application of LQPackedQs @@ -211,13 +198,13 @@ mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasCom function *(A::StridedVecOrMat, adjQ::Adjoint{<:Any,<:LQPackedQ}) Q = adjQ.parent TR = promote_type(eltype(A), eltype(Q)) - return mul!(copy_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q))) + return mul1!(copy_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q))) end function *(adjA::Adjoint{<:Any,<:StridedMatrix}, adjQ::Adjoint{<:Any,<:LQPackedQ}) A, Q = adjA.parent, adjQ.parent TR = promote_type(eltype(A), eltype(Q)) C = adjoint!(similar(A, TR, reverse(size(A))), A) - return mul!(C, adjoint(AbstractMatrix{TR}(Q))) + return mul1!(C, adjoint(AbstractMatrix{TR}(Q))) end # # (2) the inner dimension in the multiplication is the LQPackedQ's first dimension. @@ -242,7 +229,7 @@ function *(A::StridedVecOrMat, Q::LQPackedQ) else _rightappdimmismatch("columns") end - return mul!(C, AbstractMatrix{TR}(Q)) + return mul1!(C, AbstractMatrix{TR}(Q)) end function *(adjA::Adjoint{<:Any,<:StridedMatrix}, Q::LQPackedQ) A = adjA.parent @@ -255,7 +242,7 @@ function *(adjA::Adjoint{<:Any,<:StridedMatrix}, Q::LQPackedQ) else _rightappdimmismatch("rows") end - return mul!(C, AbstractMatrix{TR}(Q)) + return mul1!(C, AbstractMatrix{TR}(Q)) end _rightappdimmismatch(rowsorcols) = throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ", @@ -291,6 +278,6 @@ end function ldiv!(A::LQ{T}, B::StridedVecOrMat{T}) where T - mul!(adjoint(A.Q), ldiv!(LowerTriangular(A.L),B)) + mul2!(adjoint(A.Q), ldiv!(LowerTriangular(A.L),B)) return B end diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 3f3bf8143db3b..ae2579715ddf5 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -190,12 +190,18 @@ julia> Y mul!(C::AbstractMatrix, A::AbstractVecOrMat, B::AbstractVecOrMat) = generic_matmatmul!(C, 'N', 'N', A, B) """ - mul!(A, B) + mul1!(A, B) -Calculate the matrix-matrix product ``AB``, overwriting one of `A` or `B` (but not both), -and return the result (the overwritten argument). +Calculate the matrix-matrix product ``AB``, overwriting `A`, and return the result. """ -mul!(A, B) +mul1!(A, B) + +""" + mul2!(A, B) + +Calculate the matrix-matrix product ``AB``, overwriting `B`, and return the result. +""" +mul2!(A, B) function *(transA::Transpose{<:Any,<:AbstractMatrix}, B::AbstractMatrix) A = transA.parent diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index dae11f197942b..2fdc2473e7c80 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -329,13 +329,13 @@ function _qr(A::Union{Number,AbstractMatrix}, ::Val{false}; full::Bool = false) F = qrfact(A, Val(false)) Q, R = F.Q, F.R sQf1 = size(Q.factors, 1) - return (!full ? Array(Q) : mul!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R + return (!full ? Array(Q) : mul2!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R end function _qr(A::Union{Number, AbstractMatrix}, ::Val{true}; full::Bool = false) F = qrfact(A, Val(true)) Q, R, p = F.Q, F.R, F.p sQf1 = size(Q.factors, 1) - return (!full ? Array(Q) : mul!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R, p + return (!full ? Array(Q) : mul2!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R, p end """ @@ -506,7 +506,7 @@ AbstractMatrix{T}(Q::QRPackedQ) where {T} = QRPackedQ{T}(Q) QRCompactWYQ{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T)) AbstractMatrix{S}(Q::QRCompactWYQ{S}) where {S} = Q AbstractMatrix{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ{S}(Q) -Matrix(A::AbstractQ{T}) where {T} = mul!(A, Matrix{T}(I, size(A.factors, 1), min(size(A.factors)...))) +Matrix(A::AbstractQ{T}) where {T} = mul2!(A, Matrix{T}(I, size(A.factors, 1), min(size(A.factors)...))) Array(A::AbstractQ) = Matrix(A) size(A::Union{QR,QRCompactWY,QRPivoted}, dim::Integer) = size(getfield(A, :factors), dim) @@ -520,16 +520,16 @@ function getindex(A::AbstractQ, i::Integer, j::Integer) x[i] = 1 y = zeros(eltype(A), size(A, 2)) y[j] = 1 - return dot(x, mul!(A, y)) + return dot(x, mul2!(A, y)) end ## Multiplication by Q ### QB -mul!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = +mul2!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = LAPACK.gemqrt!('L','N',A.factors,A.T,B) -mul!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = +mul2!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = LAPACK.ormqr!('L','N',A.factors,A.τ,B) -function mul!(A::QRPackedQ, B::AbstractVecOrMat) +function mul2!(A::QRPackedQ, B::AbstractVecOrMat) mA, nA = size(A.factors) mB, nB = size(B,1), size(B,2) if mA != mB @@ -564,7 +564,7 @@ function (*)(A::AbstractQ, b::StridedVector) else throw(DimensionMismatch("vector must have length either $(size(A.factors, 1)) or $(size(A.factors, 2))")) end - mul!(Anew, bnew) + mul2!(Anew, bnew) end function (*)(A::AbstractQ, B::StridedMatrix) TAB = promote_type(eltype(A), eltype(B)) @@ -576,19 +576,19 @@ function (*)(A::AbstractQ, B::StridedMatrix) else throw(DimensionMismatch("first dimension of matrix must have size either $(size(A.factors, 1)) or $(size(A.factors, 2))")) end - mul!(Anew, Bnew) + mul2!(Anew, Bnew) end ### QcB -mul!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} = +mul2!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} = (A = adjA.parent; LAPACK.gemqrt!('L','T',A.factors,A.T,B)) -mul!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = +mul2!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = (A = adjA.parent; LAPACK.gemqrt!('L','C',A.factors,A.T,B)) -mul!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} = +mul2!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} = (A = adjA.parent; LAPACK.ormqr!('L','T',A.factors,A.τ,B)) -mul!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = +mul2!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = (A = adjA.parent; LAPACK.ormqr!('L','C',A.factors,A.τ,B)) -function mul!(adjA::Adjoint{<:Any,<:QRPackedQ}, B::AbstractVecOrMat) +function mul2!(adjA::Adjoint{<:Any,<:QRPackedQ}, B::AbstractVecOrMat) A = adjA.parent mA, nA = size(A.factors) mB, nB = size(B,1), size(B,2) @@ -616,7 +616,7 @@ end function *(adjQ::Adjoint{<:Any,<:AbstractQ}, B::StridedVecOrMat) Q = adjQ.parent TQB = promote_type(eltype(Q), eltype(B)) - return mul!(adjoint(convert(AbstractMatrix{TQB}, Q)), copy_oftype(B, TQB)) + return mul2!(adjoint(convert(AbstractMatrix{TQB}, Q)), copy_oftype(B, TQB)) end ### QBc/QcBc @@ -625,22 +625,22 @@ function *(Q::AbstractQ, adjB::Adjoint{<:Any,<:StridedVecOrMat}) TQB = promote_type(eltype(Q), eltype(B)) Bc = similar(B, TQB, (size(B, 2), size(B, 1))) adjoint!(Bc, B) - return mul!(convert(AbstractMatrix{TQB}, Q), Bc) + return mul2!(convert(AbstractMatrix{TQB}, Q), Bc) end function *(adjQ::Adjoint{<:Any,<:AbstractQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat}) Q, B = adjQ.parent, adjB.parent TQB = promote_type(eltype(Q), eltype(B)) Bc = similar(B, TQB, (size(B, 2), size(B, 1))) adjoint!(Bc, B) - return mul!(adjoint(convert(AbstractMatrix{TQB}, Q)), Bc) + return mul2!(adjoint(convert(AbstractMatrix{TQB}, Q)), Bc) end ### AQ -mul!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = +mul1!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = LAPACK.gemqrt!('R','N', B.factors, B.T, A) -mul!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = +mul1!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = LAPACK.ormqr!('R', 'N', B.factors, B.τ, A) -function mul!(A::StridedMatrix,Q::QRPackedQ) +function mul1!(A::StridedMatrix,Q::QRPackedQ) mQ, nQ = size(Q.factors) mA, nA = size(A,1), size(A,2) if nA != mQ @@ -667,19 +667,20 @@ end function (*)(A::StridedMatrix, Q::AbstractQ) TAQ = promote_type(eltype(A), eltype(Q)) - return mul!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q)) + + return mul1!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q)) end ### AQc -mul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasReal} = +mul1!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasReal} = (B = adjB.parent; LAPACK.gemqrt!('R','T',B.factors,B.T,A)) -mul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasComplex} = +mul1!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasComplex} = (B = adjB.parent; LAPACK.gemqrt!('R','C',B.factors,B.T,A)) -mul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasReal} = +mul1!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasReal} = (B = adjB.parent; LAPACK.ormqr!('R','T',B.factors,B.τ,A)) -mul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasComplex} = +mul1!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasComplex} = (B = adjB.parent; LAPACK.ormqr!('R','C',B.factors,B.τ,A)) -function mul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ}) +function mul1!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ}) Q = adjQ.parent mQ, nQ = size(Q.factors) mA, nA = size(A,1), size(A,2) @@ -711,9 +712,9 @@ function *(A::StridedMatrix, adjB::Adjoint{<:Any,<:AbstractQ}) if size(A,2) == size(B.factors, 1) AA = similar(A, TAB, size(A)) copyto!(AA, A) - return mul!(AA, adjoint(BB)) + return mul1!(AA, adjoint(BB)) elseif size(A,2) == size(B.factors,2) - return mul!([A zeros(TAB, size(A, 1), size(B.factors, 1) - size(B.factors, 2))], adjoint(BB)) + return mul1!([A zeros(TAB, size(A, 1), size(B.factors, 1) - size(B.factors, 2))], adjoint(BB)) else throw(DimensionMismatch("matrix A has dimensions $(size(A)) but matrix B has dimensions $(size(B))")) end @@ -727,20 +728,20 @@ function *(adjA::Adjoint{<:Any,<:StridedVecOrMat}, Q::AbstractQ) TAQ = promote_type(eltype(A), eltype(Q)) Ac = similar(A, TAQ, (size(A, 2), size(A, 1))) adjoint!(Ac, A) - return mul!(Ac, convert(AbstractMatrix{TAQ}, Q)) + return mul1!(Ac, convert(AbstractMatrix{TAQ}, Q)) end function *(adjA::Adjoint{<:Any,<:StridedVecOrMat}, adjQ::Adjoint{<:Any,<:AbstractQ}) A, Q = adjA.parent, adjQ.parent TAQ = promote_type(eltype(A), eltype(Q)) Ac = similar(A, TAQ, (size(A, 2), size(A, 1))) adjoint!(Ac, A) - return mul!(Ac, adjoint(convert(AbstractMatrix{TAQ}, Q))) + return mul1!(Ac, adjoint(convert(AbstractMatrix{TAQ}, Q))) end ldiv!(A::QRCompactWY{T}, b::StridedVector{T}) where {T<:BlasFloat} = - (ldiv!(UpperTriangular(A.R), view(mul!(adjoint(A.Q), b), 1:size(A, 2))); b) + (ldiv!(UpperTriangular(A.R), view(mul2!(adjoint(A.Q), b), 1:size(A, 2))); b) ldiv!(A::QRCompactWY{T}, B::StridedMatrix{T}) where {T<:BlasFloat} = - (ldiv!(UpperTriangular(A.R), view(mul!(adjoint(A.Q), B), 1:size(A, 2), 1:size(B, 2))); B) + (ldiv!(UpperTriangular(A.R), view(mul2!(adjoint(A.Q), B), 1:size(A, 2), 1:size(B, 2))); B) # Julia implementation similar to xgelsy function ldiv!(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real) where T<:BlasFloat @@ -772,7 +773,7 @@ function ldiv!(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real) where T<:BlasF rnk += 1 end C, τ = LAPACK.tzrzf!(A.factors[1:rnk,:]) - ldiv!(UpperTriangular(C[1:rnk,1:rnk]),view(mul!(adjoint(A.Q), view(B, 1:mA, 1:nrhs)), 1:rnk, 1:nrhs)) + ldiv!(UpperTriangular(C[1:rnk,1:rnk]),view(mul2!(adjoint(A.Q), view(B, 1:mA, 1:nrhs)), 1:rnk, 1:nrhs)) B[rnk+1:end,:] = zero(T) LAPACK.ormrz!('L', eltype(B)<:Complex ? 'C' : 'T', C, τ, view(B,1:nA,1:nrhs)) B[1:nA,:] = view(B, 1:nA, :)[invperm(A.p),:] @@ -786,7 +787,7 @@ function ldiv!(A::QR{T}, B::StridedMatrix{T}) where T m, n = size(A) minmn = min(m,n) mB, nB = size(B) - mul!(adjoint(A.Q), view(B, 1:m, :)) + mul2!(adjoint(A.Q), view(B, 1:m, :)) R = A.R @inbounds begin if n > m # minimum norm solution @@ -859,9 +860,7 @@ function (\)(A::Union{QR{TA},QRCompactWY{TA},QRPivoted{TA}}, B::AbstractVecOrMat X = _zeros(S, B, n) X[1:size(B, 1), :] = B - ldiv!(AA, X) - return _cut_B(X, 1:n) end diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index 2f3c32a3887b7..06ec00d34b312 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -118,8 +118,8 @@ for op in (:+, :-) end end -mul!(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = - (B = adjB.parent; mul!(full!(A), adjoint(B))) +mul1!(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = + (B = adjB.parent; mul1!(full!(A), adjoint(B))) *(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (B = adjB.parent; *(copyto!(similar(parent(A)), A), adjoint(B))) diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index 0f69a192d00c0..bb82bd18fcd34 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -457,33 +457,35 @@ fillstored!(A::UnitUpperTriangular, x) = (fillband!(A.data, x, 1, size(A,2)-1); # BlasFloat routines # ###################### -mul!(A::Tridiagonal, B::AbstractTriangular) = A*full!(B) +mul2!(A::Tridiagonal, B::AbstractTriangular) = A*full!(B) # is this necessary? + mul!(C::AbstractMatrix, A::AbstractTriangular, B::Tridiagonal) = mul!(C, copyto!(similar(parent(A)), A), B) mul!(C::AbstractMatrix, A::Tridiagonal, B::AbstractTriangular) = mul!(C, A, copyto!(similar(parent(B)), B)) mul!(C::AbstractVector, A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - (B = transB.parent; mul!(A, transpose!(C, B))) + (B = transB.parent; mul2!(A, transpose!(C, B))) mul!(C::AbstractMatrix, A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - (B = transB.parent; mul!(A, transpose!(C, B))) + (B = transB.parent; mul2!(A, transpose!(C, B))) mul!(C::AbstractMatrix, A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractVecOrMat}) = - (B = adjB.parent; mul!(A, adjoint!(C, B))) + (B = adjB.parent; mul2!(A, adjoint!(C, B))) mul!(C::AbstractVecOrMat, A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractVecOrMat}) = - (B = adjB.parent; mul!(A, adjoint!(C, B))) + (B = adjB.parent; mul2!(A, adjoint!(C, B))) + # The three methods for each op are neceesary to avoid ambiguities with definitions in matmul.jl -mul!(C::AbstractVector , A::AbstractTriangular, B::AbstractVector) = mul!(A, copyto!(C, B)) -mul!(C::AbstractMatrix , A::AbstractTriangular, B::AbstractVecOrMat) = mul!(A, copyto!(C, B)) -mul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) = mul!(A, copyto!(C, B)) +mul!(C::AbstractVector , A::AbstractTriangular, B::AbstractVector) = mul2!(A, copyto!(C, B)) +mul!(C::AbstractMatrix , A::AbstractTriangular, B::AbstractVecOrMat) = mul2!(A, copyto!(C, B)) +mul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) = mul2!(A, copyto!(C, B)) mul!(C::AbstractVector , adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVector) = - (A = adjA.parent; mul!(adjoint(A), copyto!(C, B))) + (A = adjA.parent; mul2!(adjoint(A), copyto!(C, B))) mul!(C::AbstractMatrix , adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) = - (A = adjA.parent; mul!(adjoint(A), copyto!(C, B))) + (A = adjA.parent; mul2!(adjoint(A), copyto!(C, B))) mul!(C::AbstractVecOrMat, adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) = - (A = adjA.parent; mul!(adjoint(A), copyto!(C, B))) + (A = adjA.parent; mul2!(adjoint(A), copyto!(C, B))) mul!(C::AbstractVector , transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVector) = - (A = transA.parent; mul!(transpose(A), copyto!(C, B))) + (A = transA.parent; mul2!(transpose(A), copyto!(C, B))) mul!(C::AbstractMatrix , transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) = - (A = transA.parent; mul!(transpose(A), copyto!(C, B))) + (A = transA.parent; mul2!(transpose(A), copyto!(C, B))) mul!(C::AbstractVecOrMat, transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) = - (A = transA.parent; mul!(transpose(A), copyto!(C, B))) + (A = transA.parent; mul2!(transpose(A), copyto!(C, B))) mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B)) mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B)) mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B)) @@ -491,40 +493,39 @@ mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{< mul!(C::AbstractVector, A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B))) mul!(C::AbstractVector, A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B))) - for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'), (:UnitLowerTriangular, 'L', 'U'), (:UpperTriangular, 'U', 'N'), (:UnitUpperTriangular, 'U', 'U')) @eval begin # Vector multiplication - mul!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} = + mul2!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} = BLAS.trmv!($uploc, 'N', $isunitc, A.data, b) - mul!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasFloat} = + mul2!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasFloat} = (A = transA.parent; BLAS.trmv!($uploc, 'T', $isunitc, A.data, b)) - mul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasReal} = + mul2!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasReal} = (A = adjA.parent; BLAS.trmv!($uploc, 'T', $isunitc, A.data, b)) - mul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasComplex} = + mul2!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasComplex} = (A = adjA.parent; BLAS.trmv!($uploc, 'C', $isunitc, A.data, b)) # Matrix multiplication - mul!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} = + mul2!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} = BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), A.data, B) - mul!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} = + mul1!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} = BLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), B.data, A) - mul!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasFloat} = + mul2!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasFloat} = (A = transA.parent; BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B)) - mul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasComplex} = + mul2!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasComplex} = (A = adjA.parent; BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), A.data, B)) - mul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasReal} = + mul2!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasReal} = (A = adjA.parent; BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B)) - mul!(A::StridedMatrix{T}, transB::Transpose{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasFloat} = + mul1!(A::StridedMatrix{T}, transB::Transpose{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasFloat} = (B = transB.parent; BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)) - mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasComplex} = + mul1!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasComplex} = (B = adjB.parent; BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), B.data, A)) - mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasReal} = + mul1!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasReal} = (B = adjB.parent; BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)) # Left division @@ -660,7 +661,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular), end ## Generic triangular multiplication -function mul!(A::UpperTriangular, B::StridedVecOrMat) +function mul2!(A::UpperTriangular, B::StridedVecOrMat) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) @@ -676,7 +677,8 @@ function mul!(A::UpperTriangular, B::StridedVecOrMat) end B end -function mul!(A::UnitUpperTriangular, B::StridedVecOrMat) + +function mul2!(A::UnitUpperTriangular, B::StridedVecOrMat) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) @@ -693,7 +695,7 @@ function mul!(A::UnitUpperTriangular, B::StridedVecOrMat) B end -function mul!(A::LowerTriangular, B::StridedVecOrMat) +function mul2!(A::LowerTriangular, B::StridedVecOrMat) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) @@ -709,7 +711,7 @@ function mul!(A::LowerTriangular, B::StridedVecOrMat) end B end -function mul!(A::UnitLowerTriangular, B::StridedVecOrMat) +function mul2!(A::UnitLowerTriangular, B::StridedVecOrMat) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) @@ -726,7 +728,7 @@ function mul!(A::UnitLowerTriangular, B::StridedVecOrMat) B end -function mul!(adjA::Adjoint{<:Any,<:UpperTriangular}, B::StridedVecOrMat) +function mul2!(adjA::Adjoint{<:Any,<:UpperTriangular}, B::StridedVecOrMat) A = adjA.parent m, n = size(B, 1), size(B, 2) if m != size(A, 1) @@ -743,7 +745,8 @@ function mul!(adjA::Adjoint{<:Any,<:UpperTriangular}, B::StridedVecOrMat) end B end -function mul!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat) + +function mul2!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat) A = adjA.parent m, n = size(B, 1), size(B, 2) if m != size(A, 1) @@ -761,7 +764,7 @@ function mul!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat) B end -function mul!(adjA::Adjoint{<:Any,<:LowerTriangular}, B::StridedVecOrMat) +function mul2!(adjA::Adjoint{<:Any,<:LowerTriangular}, B::StridedVecOrMat) A = adjA.parent m, n = size(B, 1), size(B, 2) if m != size(A, 1) @@ -778,7 +781,7 @@ function mul!(adjA::Adjoint{<:Any,<:LowerTriangular}, B::StridedVecOrMat) end B end -function mul!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat) +function mul2!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat) A = adjA.parent m, n = size(B, 1), size(B, 2) if m != size(A, 1) @@ -796,7 +799,7 @@ function mul!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat) B end -function mul!(transA::Transpose{<:Any,<:UpperTriangular}, B::StridedVecOrMat) +function mul2!(transA::Transpose{<:Any,<:UpperTriangular}, B::StridedVecOrMat) A = transA.parent m, n = size(B, 1), size(B, 2) if m != size(A, 1) @@ -813,7 +816,7 @@ function mul!(transA::Transpose{<:Any,<:UpperTriangular}, B::StridedVecOrMat) end B end -function mul!(transA::Transpose{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat) +function mul2!(transA::Transpose{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat) A = transA.parent m, n = size(B, 1), size(B, 2) if m != size(A, 1) @@ -831,7 +834,7 @@ function mul!(transA::Transpose{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat B end -function mul!(transA::Transpose{<:Any,<:LowerTriangular}, B::StridedVecOrMat) +function mul2!(transA::Transpose{<:Any,<:LowerTriangular}, B::StridedVecOrMat) A = transA.parent m, n = size(B, 1), size(B, 2) if m != size(A, 1) @@ -848,7 +851,7 @@ function mul!(transA::Transpose{<:Any,<:LowerTriangular}, B::StridedVecOrMat) end B end -function mul!(transA::Transpose{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat) +function mul2!(transA::Transpose{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat) A = transA.parent m, n = size(B, 1), size(B, 2) if m != size(A, 1) @@ -866,7 +869,7 @@ function mul!(transA::Transpose{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat B end -function mul!(A::StridedMatrix, B::UpperTriangular) +function mul1!(A::StridedMatrix, B::UpperTriangular) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) @@ -882,7 +885,7 @@ function mul!(A::StridedMatrix, B::UpperTriangular) end A end -function mul!(A::StridedMatrix, B::UnitUpperTriangular) +function mul1!(A::StridedMatrix, B::UnitUpperTriangular) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) @@ -899,7 +902,7 @@ function mul!(A::StridedMatrix, B::UnitUpperTriangular) A end -function mul!(A::StridedMatrix, B::LowerTriangular) +function mul1!(A::StridedMatrix, B::LowerTriangular) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) @@ -915,7 +918,7 @@ function mul!(A::StridedMatrix, B::LowerTriangular) end A end -function mul!(A::StridedMatrix, B::UnitLowerTriangular) +function mul1!(A::StridedMatrix, B::UnitLowerTriangular) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) @@ -932,7 +935,7 @@ function mul!(A::StridedMatrix, B::UnitLowerTriangular) A end -function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UpperTriangular}) +function mul1!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UpperTriangular}) B = adjB.parent m, n = size(A) if size(B, 1) != n @@ -949,7 +952,8 @@ function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UpperTriangular}) end A end -function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitUpperTriangular}) + +function mul1!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitUpperTriangular}) B = adjB.parent m, n = size(A) if size(B, 1) != n @@ -967,7 +971,7 @@ function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitUpperTriangular}) A end -function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:LowerTriangular}) +function mul1!(A::StridedMatrix, adjB::Adjoint{<:Any,<:LowerTriangular}) B = adjB.parent m, n = size(A) if size(B, 1) != n @@ -984,7 +988,8 @@ function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:LowerTriangular}) end A end -function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitLowerTriangular}) + +function mul1!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitLowerTriangular}) B = adjB.parent m, n = size(A) if size(B, 1) != n @@ -1002,7 +1007,7 @@ function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitLowerTriangular}) A end -function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UpperTriangular}) +function mul1!(A::StridedMatrix, transB::Transpose{<:Any,<:UpperTriangular}) B = transB.parent m, n = size(A) if size(B, 1) != n @@ -1019,7 +1024,7 @@ function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UpperTriangular}) end A end -function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitUpperTriangular}) +function mul1!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitUpperTriangular}) B = transB.parent m, n = size(A) if size(B, 1) != n @@ -1037,7 +1042,7 @@ function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitUpperTriangular}) A end -function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:LowerTriangular}) +function mul1!(A::StridedMatrix, transB::Transpose{<:Any,<:LowerTriangular}) B = transB.parent m, n = size(A) if size(B, 1) != n @@ -1054,7 +1059,8 @@ function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:LowerTriangular}) end A end -function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitLowerTriangular}) + +function mul1!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitLowerTriangular}) B = transB.parent m, n = size(A) if size(B, 1) != n @@ -1136,7 +1142,7 @@ function naivesub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector end # in the following transpose and conjugate transpose naive substitution variants, # accumulating in z rather than b[j] significantly improves performance as of Dec 2015 -function ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector = b) +function ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector) A = transA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1152,7 +1158,9 @@ function ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector, x: end x end -function ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector = b) +ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector) = ldiv!(transA, b, b) + +function ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector) A = transA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1167,7 +1175,9 @@ function ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector end x end -function ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector = b) +ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector) = ldiv!(transA, b, b) + +function ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector) A = transA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1183,7 +1193,9 @@ function ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector, x: end x end -function ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector = b) +ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector) = ldiv!(transA, b, b) + +function ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector) A = transA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1198,7 +1210,9 @@ function ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector end x end -function ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector = b) +ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector) = ldiv!(transA, b, b) + +function ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector) A = adjA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1214,7 +1228,9 @@ function ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::Abs end x end -function ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector = b) +ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector) = ldiv!(adjA, b, b) + +function ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector) A = adjA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1229,7 +1245,9 @@ function ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x: end x end -function ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector = b) +ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector) = ldiv!(adjA, b, b) + +function ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector) A = adjA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1245,7 +1263,9 @@ function ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector, x::Abs end x end -function ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector = b) +ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector) = ldiv!(adjA, b, b) + +function ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector) A = adjA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1260,6 +1280,7 @@ function ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x: end x end +ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector) = ldiv!(adjA, b, b) function rdiv!(A::StridedMatrix, B::UpperTriangular) m, n = size(A) @@ -1467,14 +1488,14 @@ function rdiv!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitLowerTriangular}) A end -mul!(adjA::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) = - (A = adjA.parent; UpperTriangular(mul!(adjoint(A), triu!(B.data)))) -mul!(adjA::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) = - (A = adjA.parent; LowerTriangular(mul!(adjoint(A), tril!(B.data)))) -mul!(transA::Transpose{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) = - (A = transA.parent; UpperTriangular(mul!(transpose(A), triu!(B.data)))) -mul!(transA::Transpose{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) = - (A = transA.parent; LowerTriangular(mul!(transpose(A), tril!(B.data)))) +mul2!(adjA::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) = + (A = adjA.parent; UpperTriangular(mul2!(adjoint(A), triu!(B.data)))) +mul2!(adjA::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) = + (A = adjA.parent; LowerTriangular(mul2!(adjoint(A), tril!(B.data)))) +mul2!(transA::Transpose{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) = + (A = transA.parent; UpperTriangular(mul2!(transpose(A), triu!(B.data)))) +mul2!(transA::Transpose{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) = + (A = transA.parent; LowerTriangular(mul2!(transpose(A), tril!(B.data)))) ldiv!(adjA::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) = (A = adjA.parent; UpperTriangular(ldiv!(adjoint(A), triu!(B.data)))) ldiv!(adjA::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) = @@ -1489,14 +1510,14 @@ rdiv!(A::UpperTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) = rdiv!(A::LowerTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) = LowerTriangular(rdiv!(tril!(A.data), B)) -mul!(A::UpperTriangular, adjB::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) = - (B = adjB.parent; UpperTriangular(mul!(triu!(A.data), adjoint(B)))) -mul!(A::LowerTriangular, adjB::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) = - (B = adjB.parent; LowerTriangular(mul!(tril!(A.data), adjoint(B)))) -mul!(A::UpperTriangular, transB::Transpose{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) = - (B = transB.parent; UpperTriangular(mul!(triu!(A.data), transpose(B)))) -mul!(A::LowerTriangular, transB::Transpose{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) = - (B = transB.parent; LowerTriangular(mul!(tril!(A.data), transpose(B)))) +mul1!(A::UpperTriangular, adjB::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) = + (B = adjB.parent; UpperTriangular(mul1!(triu!(A.data), adjoint(B)))) +mul1!(A::LowerTriangular, adjB::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) = + (B = adjB.parent; LowerTriangular(mul1!(tril!(A.data), adjoint(B)))) +mul1!(A::UpperTriangular, transB::Transpose{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) = + (B = transB.parent; UpperTriangular(mul1!(triu!(A.data), transpose(B)))) +mul1!(A::LowerTriangular, transB::Transpose{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) = + (B = transB.parent; LowerTriangular(mul1!(tril!(A.data), transpose(B)))) rdiv!(A::UpperTriangular, adjB::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) = (B = adjB.parent; UpperTriangular(rdiv!(triu!(A.data), adjoint(B)))) rdiv!(A::LowerTriangular, adjB::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) = @@ -1515,47 +1536,47 @@ rdiv!(A::LowerTriangular, transB::Transpose{<:Any,<:Union{UpperTriangular,UnitUp ## Some Triangular-Triangular cases. We might want to write taylored methods ## for these cases, but I'm not sure it is worth it. -(*)(A::Union{Tridiagonal,SymTridiagonal}, B::AbstractTriangular) = mul!(Matrix(A), B) +(*)(A::Union{Tridiagonal,SymTridiagonal}, B::AbstractTriangular) = mul1!(Matrix(A), B) -for (f1, f2) in ((:*, :mul!), (:\, :ldiv!)) +for (f, f2!) in ((:*, :mul2!), (:\, :ldiv!)) @eval begin - function ($f1)(A::LowerTriangular, B::LowerTriangular) - TAB = typeof(($f1)(zero(eltype(A)), zero(eltype(B))) + - ($f1)(zero(eltype(A)), zero(eltype(B)))) + function ($f)(A::LowerTriangular, B::LowerTriangular) + TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) + + ($f)(zero(eltype(A)), zero(eltype(B)))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - return LowerTriangular($f2(convert(AbstractMatrix{TAB}, A), BB)) + return LowerTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB)) end - function $(f1)(A::UnitLowerTriangular, B::LowerTriangular) + function $(f)(A::UnitLowerTriangular, B::LowerTriangular) TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) + (*)(zero(eltype(A)), zero(eltype(B)))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - return LowerTriangular($f2(convert(AbstractMatrix{TAB}, A), BB)) + return LowerTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB)) end - function ($f1)(A::UpperTriangular, B::UpperTriangular) - TAB = typeof(($f1)(zero(eltype(A)), zero(eltype(B))) + - ($f1)(zero(eltype(A)), zero(eltype(B)))) + function ($f)(A::UpperTriangular, B::UpperTriangular) + TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) + + ($f)(zero(eltype(A)), zero(eltype(B)))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - return UpperTriangular($f2(convert(AbstractMatrix{TAB}, A), BB)) + return UpperTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB)) end - function ($f1)(A::UnitUpperTriangular, B::UpperTriangular) + function ($f)(A::UnitUpperTriangular, B::UpperTriangular) TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) + (*)(zero(eltype(A)), zero(eltype(B)))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - return UpperTriangular($f2(convert(AbstractMatrix{TAB}, A), BB)) + return UpperTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB)) end end end for (ipop, op, xformtype, xformop) in ( - (:mul!, :*, :Adjoint, :adjoint), - (:mul!, :*, :Transpose, :transpose), + (:mul2!, :*, :Adjoint, :adjoint), + (:mul2!, :*, :Transpose, :transpose), (:ldiv!, :\, :Adjoint, :adjoint), (:ldiv!, :\, :Transpose, :transpose)) @eval begin @@ -1627,8 +1648,8 @@ function (/)(A::UpperTriangular, B::UnitUpperTriangular) end for (ipop, op, xformtype, xformop) in ( - (:mul!, :*, :Adjoint, :adjoint), - (:mul!, :*, :Transpose, :transpose), + (:mul1!, :*, :Adjoint, :adjoint), + (:mul1!, :*, :Transpose, :transpose), (:rdiv!, :/, :Adjoint, :adjoint), (:rdiv!, :/, :Transpose, :transpose)) @eval begin @@ -1671,26 +1692,25 @@ for (ipop, op, xformtype, xformop) in ( end ## The general promotion methods - function *(A::AbstractTriangular, B::AbstractTriangular) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - mul!(convert(AbstractArray{TAB}, A), BB) + mul2!(convert(AbstractArray{TAB}, A), BB) end function *(adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractTriangular) A = adjA.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - mul!(adjoint(convert(AbstractArray{TAB}, A)), BB) + mul2!(adjoint(convert(AbstractArray{TAB}, A)), BB) end function *(transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractTriangular) A = transA.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - mul!(transpose(convert(AbstractArray{TAB}, A)), BB) + mul2!(transpose(convert(AbstractArray{TAB}, A)), BB) end function *(A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractTriangular}) @@ -1698,14 +1718,14 @@ function *(A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractTriangular}) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) copyto!(AA, A) - mul!(AA, adjoint(convert(AbstractArray{TAB}, B))) + mul1!(AA, adjoint(convert(AbstractArray{TAB}, B))) end function *(A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractTriangular}) B = transB.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) copyto!(AA, A) - mul!(AA, transpose(convert(AbstractArray{TAB}, B))) + mul1!(AA, transpose(convert(AbstractArray{TAB}, B))) end for mat in (:AbstractVector, :AbstractMatrix) @@ -1715,21 +1735,21 @@ for mat in (:AbstractVector, :AbstractMatrix) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - mul!(convert(AbstractArray{TAB}, A), BB) + mul2!(convert(AbstractArray{TAB}, A), BB) end function *(adjA::Adjoint{<:Any,<:AbstractTriangular}, B::$mat) A = adjA.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - mul!(adjoint(convert(AbstractArray{TAB}, A)), BB) + mul2!(adjoint(convert(AbstractArray{TAB}, A)), BB) end function *(transA::Transpose{<:Any,<:AbstractTriangular}, B::$mat) A = transA.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) - mul!(transpose(convert(AbstractArray{TAB}, A)), BB) + mul2!(transpose(convert(AbstractArray{TAB}, A)), BB) end end ### Left division with triangle to the left hence rhs cannot be transposed. No quotients. @@ -1831,21 +1851,21 @@ function *(A::AbstractMatrix, B::AbstractTriangular) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) copyto!(AA, A) - mul!(AA, convert(AbstractArray{TAB}, B)) + mul1!(AA, convert(AbstractArray{TAB}, B)) end function *(A::AbstractMatrix, adjB::Adjoint{<:Any,<:AbstractTriangular}) B = adjB.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) copyto!(AA, A) - mul!(AA, adjoint(convert(AbstractArray{TAB}, B))) + mul1!(AA, adjoint(convert(AbstractArray{TAB}, B))) end function *(A::AbstractMatrix, transB::Transpose{<:Any,<:AbstractTriangular}) B = transB.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) copyto!(AA, A) - mul!(AA, transpose(convert(AbstractArray{TAB}, B))) + mul1!(AA, transpose(convert(AbstractArray{TAB}, B))) end # ambiguity resolution with definitions in linalg/rowvector.jl *(v::AdjointAbsVec, A::AbstractTriangular) = adjoint(adjoint(A) * v.parent) diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index 6c467d293f1f0..d17c25fbd4e11 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -3,7 +3,7 @@ module TestDiagonal using Test, LinearAlgebra, SparseArrays, Random -using LinearAlgebra: mul!, ldiv!, rdiv!, BlasFloat, BlasComplex, SingularException +using LinearAlgebra: mul!, mul1!, mul2!, ldiv!, rdiv!, BlasFloat, BlasComplex, SingularException n=12 #Size of matrix problem to test srand(1) @@ -147,9 +147,9 @@ srand(1) if relty <: BlasFloat b = rand(elty,n,n) b = sparse(b) - @test mul!(copy(D), copy(b)) ≈ Array(D)*Array(b) - @test mul!(transpose(copy(D)), copy(b)) ≈ transpose(Array(D))*Array(b) - @test mul!(adjoint(copy(D)), copy(b)) ≈ Array(D)'*Array(b) + @test mul2!(copy(D), copy(b)) ≈ Array(D)*Array(b) + @test mul2!(transpose(copy(D)), copy(b)) ≈ transpose(Array(D))*Array(b) + @test mul2!(adjoint(copy(D)), copy(b)) ≈ Array(D)'*Array(b) end end @@ -178,13 +178,13 @@ srand(1) VV = Array(D) DD = copy(D) r = VV * Matrix(D) - @test Array(mul!(VV, DD)) ≈ r ≈ Array(D)*Array(D) + @test Array(mul1!(VV, DD)) ≈ r ≈ Array(D)*Array(D) DD = copy(D) r = VV * transpose(Array(D)) - @test Array(mul!(VV, transpose(DD))) ≈ r + @test Array(mul1!(VV, transpose(DD))) ≈ r DD = copy(D) r = VV * Array(D)' - @test Array(mul!(VV, adjoint(DD))) ≈ r + @test Array(mul1!(VV, adjoint(DD))) ≈ r end @testset "triu/tril" begin @test istriu(D) @@ -348,9 +348,12 @@ end end let D1 = Diagonal(rand(5)), D2 = Diagonal(rand(5)) - @test_throws MethodError mul!(D1,D2) - @test_throws MethodError mul!(transpose(D1),D2) - @test_throws MethodError mul!(adjoint(D1),D2) + @test LinearAlgebra.mul1!(copy(D1),D2) == D1*D2 + @test LinearAlgebra.mul2!(D1,copy(D2)) == D1*D2 + @test LinearAlgebra.mul1!(copy(D1),transpose(D2)) == D1*transpose(D2) + @test LinearAlgebra.mul2!(transpose(D1),copy(D2)) == transpose(D1)*D2 + @test LinearAlgebra.mul1!(copy(D1),adjoint(D2)) == D1*adjoint(D2) + @test LinearAlgebra.mul2!(adjoint(D1),copy(D2)) == adjoint(D1)*D2 end @testset "multiplication of QR Q-factor and Diagonal (#16615 spot test)" begin @@ -358,7 +361,7 @@ end Q = qrfact(randn(5, 5)).Q @test D * Q' == Array(D) * Q' Q = qrfact(randn(5, 5), Val(true)).Q - @test_throws MethodError mul!(Q, D) + @test_throws ArgumentError mul2!(Q, D) end @testset "block diagonal matrices" begin diff --git a/stdlib/LinearAlgebra/test/givens.jl b/stdlib/LinearAlgebra/test/givens.jl index fd722d8aa4868..df15b0a088926 100644 --- a/stdlib/LinearAlgebra/test/givens.jl +++ b/stdlib/LinearAlgebra/test/givens.jl @@ -3,6 +3,7 @@ module TestGivens using Test, LinearAlgebra, Random +using LinearAlgebra: mul1!, mul2! # Test givens rotations @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) @@ -17,11 +18,11 @@ using Test, LinearAlgebra, Random for j = 1:8 for i = j+2:10 G, _ = givens(A, j+1, i, j) - mul!(G, A) - mul!(A, adjoint(G)) - mul!(G, R) + mul2!(G, A) + mul1!(A, adjoint(G)) + mul2!(G, R) - @test mul!(G,Matrix{elty}(I, 10, 10)) == [G[i,j] for i=1:10,j=1:10] + @test mul2!(G,Matrix{elty}(I, 10, 10)) == [G[i,j] for i=1:10,j=1:10] @testset "transposes" begin @test copy(G')*G*Matrix(elty(1)I, 10, 10) ≈ Matrix(I, 10, 10) @@ -34,8 +35,8 @@ using Test, LinearAlgebra, Random @test_throws ArgumentError givens(A, 3, 3, 2) @test_throws ArgumentError givens(one(elty),zero(elty),2,2) G, _ = givens(one(elty),zero(elty),11,12) - @test_throws DimensionMismatch mul!(G, A) - @test_throws DimensionMismatch mul!(A, adjoint(G)) + @test_throws DimensionMismatch mul2!(G, A) + @test_throws DimensionMismatch mul1!(A, adjoint(G)) @test abs.(A) ≈ abs.(hessfact(Ac).H) @test norm(R*Matrix{elty}(I, 10, 10)) ≈ one(elty) diff --git a/stdlib/LinearAlgebra/test/lq.jl b/stdlib/LinearAlgebra/test/lq.jl index dc63e77a0441c..8b51b496f09f9 100644 --- a/stdlib/LinearAlgebra/test/lq.jl +++ b/stdlib/LinearAlgebra/test/lq.jl @@ -3,7 +3,7 @@ module TestLQ using Test, LinearAlgebra, Random -using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, mul! +using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, mul1!, mul2! n = 10 @@ -21,7 +21,7 @@ breal = randn(n,2)/2 bimg = randn(n,2)/2 # helper functions to unambiguously recover explicit forms of an LQPackedQ -squareQ(Q::LinearAlgebra.LQPackedQ) = (n = size(Q.factors, 2); mul!(Q, Matrix{eltype(Q)}(I, n, n))) +squareQ(Q::LinearAlgebra.LQPackedQ) = (n = size(Q.factors, 2); mul2!(Q, Matrix{eltype(Q)}(I, n, n))) rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q) @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64) @@ -54,8 +54,6 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q) @test size(lqa.Q,3) == 1 @test_throws ErrorException lqa.Z @test Array(copy(adjoint(lqa))) ≈ a' - @test lqa * lqa' ≈ a * a' - @test lqa' * lqa ≈ a' * a @test q*squareQ(q)' ≈ Matrix(I, n, n) @test l*q ≈ a @test Array(lqa) ≈ a @@ -99,9 +97,9 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q) l,q = lqa.L, lqa.Q @test rectangularQ(q)*rectangularQ(q)' ≈ Matrix(I, n1, n1) @test squareQ(q)'*squareQ(q) ≈ Matrix(I, n1, n1) - @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1),q) - @test mul!(adjoint(q), rectangularQ(q)) ≈ Matrix(I, n1, n1) - @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1), adjoint(q)) + @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1),q) + @test mul2!(adjoint(q), rectangularQ(q)) ≈ Matrix(I, n1, n1) + @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1), adjoint(q)) @test_throws BoundsError size(q,-1) end end @@ -149,7 +147,7 @@ end function getqs(F::LinearAlgebra.LQ) implicitQ = F.Q sq = size(implicitQ.factors, 2) - explicitQ = mul!(implicitQ, Matrix{eltype(implicitQ)}(I, sq, sq)) + explicitQ = mul2!(implicitQ, Matrix{eltype(implicitQ)}(I, sq, sq)) return implicitQ, explicitQ end @@ -190,7 +188,7 @@ end @testset "postmultiplication with / right-application of LQPackedQ (#23779)" begin function getqs(F::LinearAlgebra.LQ) implicitQ = F.Q - explicitQ = mul!(implicitQ, Matrix{eltype(implicitQ)}(I, size(implicitQ)...)) + explicitQ = mul2!(implicitQ, Matrix{eltype(implicitQ)}(I, size(implicitQ)...)) return implicitQ, explicitQ end # for any shape m-by-n of LQ-factored matrix, where Q is an LQPackedQ diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl index b3ace0f81089c..99385030b556d 100644 --- a/stdlib/LinearAlgebra/test/qr.jl +++ b/stdlib/LinearAlgebra/test/qr.jl @@ -3,7 +3,7 @@ module TestQR using Test, LinearAlgebra, Random -using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, mul! +using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, mul1!, mul2! n = 10 @@ -21,7 +21,7 @@ breal = randn(n,2)/2 bimg = randn(n,2)/2 # helper functions to unambiguously recover explicit forms of an implicit QR Q -squareQ(Q::LinearAlgebra.AbstractQ) = (sq = size(Q.factors, 1); mul!(Q, Matrix{eltype(Q)}(I, sq, sq))) +squareQ(Q::LinearAlgebra.AbstractQ) = (sq = size(Q.factors, 1); mul2!(Q, Matrix{eltype(Q)}(I, sq, sq))) rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) @@ -136,20 +136,20 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) a = raw_a qrpa = factorize(a[:,1:n1]) q, r = qrpa.Q, qrpa.R - @test mul!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n) - @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1),q) - @test mul!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n) - @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1), adjoint(q)) + @test mul1!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n) + @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1),q) + @test mul1!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n) + @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1), adjoint(q)) @test_throws BoundsError size(q,-1) - @test_throws DimensionMismatch LinearAlgebra.mul!(q,zeros(eltya,n1+1)) - @test_throws DimensionMismatch LinearAlgebra.mul!(adjoint(q), zeros(eltya,n1+1)) + @test_throws DimensionMismatch LinearAlgebra.mul2!(q,zeros(eltya,n1+1)) + @test_throws DimensionMismatch LinearAlgebra.mul2!(adjoint(q), zeros(eltya,n1+1)) qra = qrfact(a[:,1:n1], Val(false)) q, r = qra.Q, qra.R - @test mul!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n) - @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1),q) - @test mul!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n) - @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1),adjoint(q)) + @test mul1!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n) + @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1),q) + @test mul1!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n) + @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1),adjoint(q)) @test_throws BoundsError size(q,-1) @test_throws DimensionMismatch q * Matrix{Int8}(I, n+4, n+4) end diff --git a/stdlib/LinearAlgebra/test/special.jl b/stdlib/LinearAlgebra/test/special.jl index 8799c298f5e3e..9563ce13507ca 100644 --- a/stdlib/LinearAlgebra/test/special.jl +++ b/stdlib/LinearAlgebra/test/special.jl @@ -3,7 +3,7 @@ module TestSpecial using Test, LinearAlgebra, SparseArrays, Random -using LinearAlgebra: mul! +using LinearAlgebra: mul1! n= 10 #Size of matrix to test srand(1) @@ -118,10 +118,10 @@ end b = rand(n,n) qrb = qrfact(b,Val(true)) @test *(atri, adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q' - @test mul!(copy(atri), adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q' + @test mul1!(copy(atri), adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q' qrb = qrfact(b,Val(false)) @test *(atri, adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q' - @test mul!(copy(atri), adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q' + @test mul1!(copy(atri), adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q' end end diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl index 08469d27fd589..c2e6420af3734 100644 --- a/stdlib/LinearAlgebra/test/triangular.jl +++ b/stdlib/LinearAlgebra/test/triangular.jl @@ -6,7 +6,7 @@ debug = false using Test, LinearAlgebra, SparseArrays, Random using LinearAlgebra: BlasFloat, errorbounds, full!, naivesub!, transpose!, UnitUpperTriangular, UnitLowerTriangular, - mul!, rdiv! + mul!, rdiv!, mul1!, mul2! debug && println("Triangular matrices") @@ -322,7 +322,7 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo if !(eltyB in (BigFloat, Complex{BigFloat})) # rand does not support BigFloat and Complex{BigFloat} as of Dec 2015 Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1)) - @test mul!(Tri,copy(A1)) ≈ Tri*Matrix(A1) + @test mul2!(Tri,copy(A1)) ≈ Tri*Matrix(A1) end # Triangular-dense Matrix/vector multiplication @@ -360,12 +360,12 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo end #error handling Ann, Bmm, bm = A1, Matrix{eltyB}(uninitialized, n+1, n+1), Vector{eltyB}(uninitialized, n+1) - @test_throws DimensionMismatch mul!(Ann, bm) - @test_throws DimensionMismatch mul!(Bmm, Ann) - @test_throws DimensionMismatch mul!(transpose(Ann), bm) - @test_throws DimensionMismatch mul!(adjoint(Ann), bm) - @test_throws DimensionMismatch mul!(Bmm, adjoint(Ann)) - @test_throws DimensionMismatch mul!(Bmm, transpose(Ann)) + @test_throws DimensionMismatch mul2!(Ann, bm) + @test_throws DimensionMismatch mul1!(Bmm, Ann) + @test_throws DimensionMismatch mul2!(transpose(Ann), bm) + @test_throws DimensionMismatch mul2!(adjoint(Ann), bm) + @test_throws DimensionMismatch mul1!(Bmm, adjoint(Ann)) + @test_throws DimensionMismatch mul1!(Bmm, transpose(Ann)) # ... and division @test A1\B[:,1] ≈ Matrix(A1)\B[:,1] diff --git a/stdlib/SparseArrays/src/deprecated.jl b/stdlib/SparseArrays/src/deprecated.jl index 829d6c9401a3d..87be3f449529c 100644 --- a/stdlib/SparseArrays/src/deprecated.jl +++ b/stdlib/SparseArrays/src/deprecated.jl @@ -222,7 +222,6 @@ end import Base: asyncmap @deprecate asyncmap(f, s::AbstractSparseArray...; kwargs...) sparse(asyncmap(f, map(Array, s)...; kwargs...)) - # END 0.7 deprecations # BEGIN 1.0 deprecations diff --git a/stdlib/SuiteSparse/src/deprecated.jl b/stdlib/SuiteSparse/src/deprecated.jl index 2c36db610d028..14bdd94e7c2ac 100644 --- a/stdlib/SuiteSparse/src/deprecated.jl +++ b/stdlib/SuiteSparse/src/deprecated.jl @@ -30,12 +30,12 @@ end LinearAlgebra.ldiv!(X, transpose(lu), B) LinearAlgebra.Ac_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} = LinearAlgebra.ldiv!(X, adjoint(lu), B) - LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(lu, B) - LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(transpose(lu), B) - LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(adjoint(lu), B) - LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(lu, B) - LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(transpose(lu), B) - LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(adjoint(lu), B) + LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, lu, copy(B)) + LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, transpose(lu), copy(B)) + LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, adjoint(lu), copy(B)) + LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, lu, copy(B)) + LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, transpose(lu), copy(B)) + LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, adjoint(lu), copy(B)) end # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from src/spqr.jl, to deprecate diff --git a/stdlib/SuiteSparse/src/spqr.jl b/stdlib/SuiteSparse/src/spqr.jl index 742342d7efa34..4245c264d8da6 100644 --- a/stdlib/SuiteSparse/src/spqr.jl +++ b/stdlib/SuiteSparse/src/spqr.jl @@ -198,7 +198,7 @@ LinearAlgebra.qrfact(A::SparseMatrixCSC; tol = _default_tol(A)) = qrfact(A, Val{ LinearAlgebra.qr(A::SparseMatrixCSC; tol = _default_tol(A)) = qr(A, Val{true}, tol = tol) -function LinearAlgebra.mul!(Q::QRSparseQ, A::StridedVecOrMat) +function LinearAlgebra.mul2!(Q::QRSparseQ, A::StridedVecOrMat) if size(A, 1) != size(Q, 1) throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) end @@ -213,7 +213,7 @@ function LinearAlgebra.mul!(Q::QRSparseQ, A::StridedVecOrMat) return A end -function LinearAlgebra.mul!(A::StridedMatrix, Q::QRSparseQ) +function LinearAlgebra.mul1!(A::StridedMatrix, Q::QRSparseQ) if size(A, 2) != size(Q, 1) throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) end @@ -227,7 +227,7 @@ function LinearAlgebra.mul!(A::StridedMatrix, Q::QRSparseQ) return A end -function LinearAlgebra.mul!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMat) +function LinearAlgebra.mul2!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMat) Q = adjQ.parent if size(A, 1) != size(Q, 1) throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) @@ -243,7 +243,7 @@ function LinearAlgebra.mul!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMat return A end -function LinearAlgebra.mul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRSparseQ}) +function LinearAlgebra.mul1!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRSparseQ}) Q = adjQ.parent if size(A, 2) != size(Q, 1) throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) @@ -375,7 +375,7 @@ function _ldiv_basic(F::QRSparse, B::StridedVecOrMat) X0 = view(X, 1:size(B, 1), :) # Apply Q' to B - LinearAlgebra.mul!(adjoint(F.Q), X0) + LinearAlgebra.mul2!(adjoint(F.Q), X0) # Zero out to get basic solution X[rnk + 1:end, :] = 0 diff --git a/stdlib/SuiteSparse/src/umfpack.jl b/stdlib/SuiteSparse/src/umfpack.jl index 9cc32f2261e23..c775c5c94aafa 100644 --- a/stdlib/SuiteSparse/src/umfpack.jl +++ b/stdlib/SuiteSparse/src/umfpack.jl @@ -388,6 +388,9 @@ function nnz(lu::UmfpackLU) end ### Solve with Factorization + +import LinearAlgebra.ldiv! + ldiv!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = ldiv!(B, lu, copy(B)) ldiv!(translu::Transpose{T,<:UmfpackLU{T}}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = diff --git a/stdlib/SuiteSparse/test/spqr.jl b/stdlib/SuiteSparse/test/spqr.jl index e59e1ff241026..c36545e173cd2 100644 --- a/stdlib/SuiteSparse/test/spqr.jl +++ b/stdlib/SuiteSparse/test/spqr.jl @@ -2,7 +2,7 @@ using SuiteSparse.SPQR using SuiteSparse.CHOLMOD -using LinearAlgebra: mul!, Adjoint, Transpose +using LinearAlgebra: mul1!, mul2!, Adjoint, Transpose @testset "Sparse QR" begin m, n = 100, 10 @@ -43,10 +43,10 @@ nn = 100 @test norm(R0[n + 1:end, :], 1) < 1e-12 offsizeA = Matrix{Float64}(I, m+1, m+1) - @test_throws DimensionMismatch mul!(Q, offsizeA) - @test_throws DimensionMismatch mul!(adjoint(Q), offsizeA) - @test_throws DimensionMismatch mul!(offsizeA, Q) - @test_throws DimensionMismatch mul!(offsizeA, adjoint(Q)) + @test_throws DimensionMismatch mul2!(Q, offsizeA) + @test_throws DimensionMismatch mul2!(adjoint(Q), offsizeA) + @test_throws DimensionMismatch mul1!(offsizeA, Q) + @test_throws DimensionMismatch mul1!(offsizeA, adjoint(Q)) end @testset "element type of B: $eltyB" for eltyB in (Int, Float64, Complex{Float64}) diff --git a/stdlib/SuiteSparse/test/umfpack.jl b/stdlib/SuiteSparse/test/umfpack.jl index 1e48075937895..a1be43d13dfc2 100644 --- a/stdlib/SuiteSparse/test/umfpack.jl +++ b/stdlib/SuiteSparse/test/umfpack.jl @@ -32,7 +32,7 @@ @test A*x ≈ b z = complex.(b) - x = SuiteSparse.ldiv!(lua, z) + x = LinearAlgebra.ldiv!(lua, z) @test x ≈ float([1:5;]) @test z === x y = similar(z) @@ -47,11 +47,11 @@ @test A'*x ≈ b z = complex.(b) - x = SuiteSparse.ldiv!(adjoint(lua), z) + x = LinearAlgebra.ldiv!(adjoint(lua), z) @test x ≈ float([1:5;]) @test x === z y = similar(x) - SuiteSparse.ldiv!(y, adjoint(lua), complex.(b)) + LinearAlgebra.ldiv!(y, adjoint(lua), complex.(b)) @test y ≈ x @test A'*x ≈ b @@ -59,10 +59,10 @@ @test x ≈ float([1:5;]) @test transpose(A) * x ≈ b - x = SuiteSparse.ldiv!(transpose(lua), complex.(b)) + x = LinearAlgebra.ldiv!(transpose(lua), complex.(b)) @test x ≈ float([1:5;]) y = similar(x) - SuiteSparse.ldiv!(y, transpose(lua), complex.(b)) + LinearAlgebra.ldiv!(y, transpose(lua), complex.(b)) @test y ≈ x @test transpose(A) * x ≈ b