From fed0945ff0e44e58495b281cbd30b5c7b689c955 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Tue, 21 Nov 2017 23:31:02 -0800 Subject: [PATCH] Clarify which args are modified in mul! --- base/deprecated.jl | 278 +++++++++++++++++--------- base/hashing2.jl | 2 +- base/linalg/diagonal.jl | 63 +++--- base/linalg/givens.jl | 33 +-- base/linalg/hessenberg.jl | 19 +- base/linalg/linalg.jl | 27 ++- base/linalg/lq.jl | 57 ++---- base/linalg/matmul.jl | 14 +- base/linalg/qr.jl | 71 ++++--- base/linalg/special.jl | 4 +- base/linalg/triangular.jl | 230 +++++++++++---------- stdlib/SparseArrays/test/sparse.jl | 2 +- stdlib/SuiteSparse/src/SuiteSparse.jl | 2 +- stdlib/SuiteSparse/src/deprecated.jl | 12 +- stdlib/SuiteSparse/src/spqr.jl | 12 +- stdlib/SuiteSparse/src/umfpack.jl | 13 +- stdlib/SuiteSparse/test/spqr.jl | 10 +- stdlib/SuiteSparse/test/umfpack.jl | 10 +- test/linalg/diagonal.jl | 27 +-- test/linalg/givens.jl | 14 +- test/linalg/lq.jl | 16 +- test/linalg/qr.jl | 24 +-- test/linalg/special.jl | 14 +- test/linalg/triangular.jl | 16 +- 24 files changed, 524 insertions(+), 446 deletions(-) diff --git a/base/deprecated.jl b/base/deprecated.jl index 93e97f884d544..bc4d03d3fa985 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -1255,6 +1255,7 @@ import .LinAlg: diff finalizer(f::Ptr{Cvoid}, o::Ptr{Cvoid}) = invoke(finalizer, Tuple{Ptr{Cvoid}, Any}, f, o) finalizer(f::Ptr{Cvoid}, o::Function) = invoke(finalizer, Tuple{Ptr{Cvoid}, Any}, f, o) + # Broadcast extension API (#23939) @eval Broadcast begin Base.@deprecate_binding containertype combine_styles false @@ -1708,11 +1709,11 @@ end # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/diagonal.jl, to deprecate @eval Base.LinAlg begin - @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) @@ -1728,16 +1729,19 @@ end @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) @@ -1763,7 +1767,7 @@ end # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/special.jl, to deprecate @eval Base.LinAlg begin - @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)) end @@ -1796,10 +1800,10 @@ end # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/hessenberg.jl, to deprecate @eval Base.LinAlg begin - @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)) end @@ -1859,18 +1863,18 @@ end # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/lq.jl, to deprecate @eval Base.LinAlg begin - @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) @@ -1879,25 +1883,25 @@ end # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/qr.jl, to deprecate @eval Base.LinAlg begin - @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) @@ -1994,20 +1998,20 @@ end @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) @@ -2032,30 +2036,30 @@ end @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) @@ -2065,7 +2069,7 @@ end @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)) @@ -2120,22 +2124,22 @@ for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'), (:UnitUpperTriangular, 'U', 'U')) @eval Base.LinAlg 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} mul1!(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) @@ -2196,15 +2200,92 @@ end # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/givens.jl, to deprecate @eval Base.LinAlg begin - @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)) end +# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/sparse/linalg.jl, to deprecate +@eval Base.SparseArrays begin + import Base: A_mul_B!, Ac_mul_B, Ac_mul_B!, At_mul_B, At_mul_B! + import Base: A_mul_Bc, A_mul_Bt, Ac_mul_Bc, At_mul_Bt + import Base: At_ldiv_B, Ac_ldiv_B, A_ldiv_B! + import Base.LinAlg: At_ldiv_B!, Ac_ldiv_B!, A_rdiv_B!, A_rdiv_Bc! + + using Base.LinAlg: Adjoint, Transpose + @deprecate Ac_ldiv_B(A::SparseMatrixCSC, B::RowVector) (\)(adjoint(A), B) + @deprecate At_ldiv_B(A::SparseMatrixCSC, B::RowVector) (\)(transpose(A), B) + @deprecate Ac_ldiv_B(A::SparseMatrixCSC, B::AbstractVecOrMat) (\)(adjoint(A), B) + @deprecate At_ldiv_B(A::SparseMatrixCSC, B::AbstractVecOrMat) (\)(transpose(A), B) + @deprecate A_rdiv_Bc!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where {T} rdiv!(A, adjoint(D)) + @deprecate A_rdiv_Bt!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where {T} rdiv!(A, transpose(D)) + @deprecate A_rdiv_B!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where {T} rdiv!(A, D) + @deprecate A_ldiv_B!(L::LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T} ldiv!(L, B) + @deprecate A_ldiv_B!(U::UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T} ldiv!(U, B) + @deprecate A_mul_Bt(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} (*)(A, transpose(B)) + @deprecate A_mul_Bc(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} (*)(A, adjoint(B)) + @deprecate At_mul_B(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} (*)(transpose(A), B) + @deprecate Ac_mul_B(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} (*)(adjoint(A), B) + @deprecate At_mul_Bt(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} (*)(transpose(A), transpose(B)) + @deprecate Ac_mul_Bc(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} (*)(adjoint(A), adjoint(B)) + @deprecate A_mul_B!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat) mul!(C, A, B) + @deprecate Ac_mul_B!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat) mul!(C, adjoint(A), B) + @deprecate At_mul_B!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat) mul!(C, transpose(A), B) + @deprecate A_mul_B!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat) mul!(α, A, B, β, C) + @deprecate A_mul_B(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx} (*)(A, x) + @deprecate A_mul_B(A::SparseMatrixCSC{TA,S}, B::StridedMatrix{Tx}) where {TA,S,Tx} (*)(A, B) + @deprecate Ac_mul_B!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat) mul!(α, adjoint(A), B, β, C) + @deprecate Ac_mul_B(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx} (*)(adjoint(A), x) + @deprecate Ac_mul_B(A::SparseMatrixCSC{TA,S}, B::StridedMatrix{Tx}) where {TA,S,Tx} (*)(adjoint(A), B) + @deprecate At_mul_B!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat) mul!(α, transpose(A), B, β, C) + @deprecate At_mul_B(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx} (*)(transpose(A), x) + @deprecate At_mul_B(A::SparseMatrixCSC{TA,S}, B::StridedMatrix{Tx}) where {TA,S,Tx} (*)(transpose(A), B) + @deprecate A_mul_Bt(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} (*)(A, transpose(B)) + @deprecate A_mul_Bc(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} (*)(A, adjoint(B)) + @deprecate At_mul_B(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} (*)(transpose(A), B) + @deprecate Ac_mul_B(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} (*)(adjoint(A),B) + @deprecate At_mul_Bt(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} (*)(transpose(A), transpose(B)) + @deprecate Ac_mul_Bc(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} (*)(adjoint(A), adjoint(B)) +end + +# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/sparse/sparsevector.jl, to deprecate +for isunittri in (true, false), islowertri in (true, false) + unitstr = isunittri ? "Unit" : "" + halfstr = islowertri ? "Lower" : "Upper" + tritype = :(Base.LinAlg.$(Symbol(unitstr, halfstr, "Triangular"))) + @eval Base.SparseArrays begin + using Base.LinAlg: Adjoint, Transpose + @deprecate At_ldiv_B(A::$tritype{TA,<:AbstractMatrix}, b::SparseVector{Tb}) where {TA<:Number,Tb<:Number} (\)(transpose(A), b) + @deprecate At_ldiv_B(A::$tritype{TA,<:StridedMatrix}, b::SparseVector{Tb}) where {TA<:Number,Tb<:Number} (\)(transpose(A), b) + @deprecate At_ldiv_B(A::$tritype, b::SparseVector) (\)(transpose(A), b) + @deprecate Ac_ldiv_B(A::$tritype{TA,<:AbstractMatrix}, b::SparseVector{Tb}) where {TA<:Number,Tb<:Number} (\)(adjoint(A), b) + @deprecate Ac_ldiv_B(A::$tritype{TA,<:StridedMatrix}, b::SparseVector{Tb}) where {TA<:Number,Tb<:Number} (\)(adjoint(A), b) + @deprecate Ac_ldiv_B(A::$tritype, b::SparseVector) (\)(adjoint(A), b) + @deprecate A_ldiv_B!(A::$tritype{<:Any,<:StridedMatrix}, b::SparseVector) ldiv!(A, b) + @deprecate At_ldiv_B!(A::$tritype{<:Any,<:StridedMatrix}, b::SparseVector) ldiv!(transpose(A), b) + @deprecate Ac_ldiv_B!(A::$tritype{<:Any,<:StridedMatrix}, b::SparseVector) ldiv!(adjoint(A), b) + end +end +@eval Base.SparseArrays begin + using Base.LinAlg: Adjoint, Transpose + @deprecate Ac_mul_B(A::SparseMatrixCSC, x::AbstractSparseVector) (*)(adjoint(A), x) + @deprecate At_mul_B(A::SparseMatrixCSC, x::AbstractSparseVector) (*)(transpose(A), x) + @deprecate Ac_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::StridedVector) mul!(α, adjoint(A), x, β, y) + @deprecate Ac_mul_B!(y::StridedVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty} mul!(y, adjoint(A), x) + @deprecate At_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::StridedVector) mul!(α, transpose(A), x, β, y) + @deprecate At_mul_B!(y::StridedVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty} mul!(y, transpose(A), x) + @deprecate A_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::StridedVector) mul!(α, A, x, β, y) + @deprecate A_mul_B!(y::StridedVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty} mul!(y, A, x) + @deprecate At_mul_B!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number, y::StridedVector) mul!(α, transpose(A), x, β, y) + @deprecate At_mul_B!(y::StridedVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where {Tx,Ty} mul!(y, transpose(A), x) + @deprecate At_mul_B(A::StridedMatrix{Ta}, x::AbstractSparseVector{Tx}) where {Ta,Tx} (*)(transpose(A), x) + @deprecate A_mul_B!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number, y::StridedVector) mul!(α, A, x, β, y) + @deprecate A_mul_B!(y::StridedVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where {Tx,Ty} mul!(y, A, x) +end # methods involving RowVector from base/linalg/bidiag.jl, to deprecate @eval Base.LinAlg begin @@ -2717,6 +2798,11 @@ end @deprecate substrides(s, parent, dim, I::Tuple) substrides(parent, strides(parent), I) +@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) + @deprecate lexcmp(x::AbstractArray, y::AbstractArray) cmp(x, y) @deprecate lexcmp(x::Real, y::Real) cmp(isless, x, y) @deprecate lexcmp(x::Complex, y::Complex) cmp((real(x),imag(x)), (real(y),imag(y))) 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/linalg/diagonal.jl b/base/linalg/diagonal.jl index 8e3f06e1000d1..8d0c8a67fb71c 100644 --- a/base/linalg/diagonal.jl +++ b/base/linalg/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,21 +233,9 @@ 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))) +mul1!(A::Diagonal,B::Diagonal) = Diagonal(A.diag .*= B.diag) +mul2!(A::Diagonal,B::Diagonal) = Diagonal(B.diag .= A.diag .* B.diag) + 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)) @@ -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/base/linalg/givens.jl b/base/linalg/givens.jl index f305760e8159f..5b5648a30838c 100644 --- a/base/linalg/givens.jl +++ b/base/linalg/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 """ LinAlg.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/base/linalg/hessenberg.jl b/base/linalg/hessenberg.jl index 5b4ea46f5fedc..8ba45b75a10a3 100644 --- a/base/linalg/hessenberg.jl +++ b/base/linalg/hessenberg.jl @@ -71,7 +71,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 @@ -82,31 +82,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/base/linalg/linalg.jl b/base/linalg/linalg.jl index 99d308b764426..3e774211bdd90 100644 --- a/base/linalg/linalg.jl +++ b/base/linalg/linalg.jl @@ -238,11 +238,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)). @@ -254,11 +252,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)). @@ -267,7 +278,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/base/linalg/lq.jl b/base/linalg/lq.jl index 77c3d7d11e232..685f78d52c6c0 100644 --- a/base/linalg/lq.jl +++ b/base/linalg/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.τ)) @@ -86,7 +86,7 @@ function getproperty(F::LQ, d::Symbol) end 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:") @@ -120,47 +120,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 @@ -172,14 +159,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 @@ -187,11 +174,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 @@ -209,13 +196,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. @@ -240,7 +227,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 @@ -253,7 +240,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 ", @@ -289,6 +276,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/base/linalg/matmul.jl b/base/linalg/matmul.jl index 229667b4bf451..7434ee066fd2a 100644 --- a/base/linalg/matmul.jl +++ b/base/linalg/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/base/linalg/qr.jl b/base/linalg/qr.jl index 376bdefba6940..ecebffef3290e 100644 --- a/base/linalg/qr.jl +++ b/base/linalg/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 """ @@ -504,7 +504,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) @@ -518,16 +518,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 @@ -562,7 +562,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)) @@ -574,19 +574,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) @@ -614,7 +614,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 @@ -623,22 +623,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 @@ -665,19 +665,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) @@ -709,9 +710,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 @@ -725,20 +726,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 @@ -770,7 +771,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),:] @@ -784,7 +785,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 @@ -857,9 +858,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/base/linalg/special.jl b/base/linalg/special.jl index 2f3c32a3887b7..06ec00d34b312 100644 --- a/base/linalg/special.jl +++ b/base/linalg/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/base/linalg/triangular.jl b/base/linalg/triangular.jl index bfb25b5f14e9e..69f18d57e25ab 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/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/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index bdbd652e3877e..e841b71515ae1 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -1,6 +1,6 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license -using Base.LinAlg: mul!, ldiv!, rdiv! +using Base.LinAlg: mul!, rdiv! using Base.Printf: @printf using Random diff --git a/stdlib/SuiteSparse/src/SuiteSparse.jl b/stdlib/SuiteSparse/src/SuiteSparse.jl index 6b7558310a1d0..a4b4ad8658f5d 100644 --- a/stdlib/SuiteSparse/src/SuiteSparse.jl +++ b/stdlib/SuiteSparse/src/SuiteSparse.jl @@ -5,7 +5,7 @@ __precompile__(true) module SuiteSparse import Base: \ -import Base.LinAlg: ldiv!, rdiv! +import Base.LinAlg: ldiv! ## Functions to switch to 0-based indexing to call external sparse solvers diff --git a/stdlib/SuiteSparse/src/deprecated.jl b/stdlib/SuiteSparse/src/deprecated.jl index 2bf2d022684a2..7028d5dfc12f0 100644 --- a/stdlib/SuiteSparse/src/deprecated.jl +++ b/stdlib/SuiteSparse/src/deprecated.jl @@ -30,12 +30,12 @@ end Base.LinAlg.ldiv!(X, transpose(lu), B) Base.Ac_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} = Base.LinAlg.ldiv!(X, adjoint(lu), B) - Base.A_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(lu, B) - Base.At_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(transpose(lu), B) - Base.Ac_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(adjoint(lu), B) - Base.A_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.ldiv!(lu, B) - Base.At_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.ldiv!(transpose(lu), B) - Base.Ac_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.ldiv!(adjoint(lu), B) + Base.A_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(B, lu, copy(B)) + Base.At_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(B, transpose(lu), copy(B)) + Base.Ac_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(B, adjoint(lu), copy(B)) + Base.A_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.ldiv!(B, lu, copy(B)) + Base.At_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.ldiv!(B, transpose(lu), copy(B)) + Base.Ac_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.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 acd3d37bbd2c7..a2f828659276e 100644 --- a/stdlib/SuiteSparse/src/spqr.jl +++ b/stdlib/SuiteSparse/src/spqr.jl @@ -197,7 +197,7 @@ Base.LinAlg.qrfact(A::SparseMatrixCSC; tol = _default_tol(A)) = qrfact(A, Val{tr Base.LinAlg.qr(A::SparseMatrixCSC; tol = _default_tol(A)) = qr(A, Val{true}, tol = tol) -function Base.LinAlg.mul!(Q::QRSparseQ, A::StridedVecOrMat) +function Base.LinAlg.mul2!(Q::QRSparseQ, A::StridedVecOrMat) if size(A, 1) != size(Q, 1) throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) end @@ -212,7 +212,7 @@ function Base.LinAlg.mul!(Q::QRSparseQ, A::StridedVecOrMat) return A end -function Base.LinAlg.mul!(A::StridedMatrix, Q::QRSparseQ) +function Base.LinAlg.mul1!(A::StridedMatrix, Q::QRSparseQ) if size(A, 2) != size(Q, 1) throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) end @@ -226,7 +226,7 @@ function Base.LinAlg.mul!(A::StridedMatrix, Q::QRSparseQ) return A end -function Base.LinAlg.mul!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMat) +function Base.LinAlg.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))")) @@ -242,7 +242,7 @@ function Base.LinAlg.mul!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMat) return A end -function Base.LinAlg.mul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRSparseQ}) +function Base.LinAlg.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))")) @@ -374,8 +374,8 @@ function _ldiv_basic(F::QRSparse, B::StridedVecOrMat) X0 = view(X, 1:size(B, 1), :) # Apply Q' to B - Base.LinAlg.mul!(adjoint(F.Q), X0) - + Base.LinAlg.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 371f1dc5d91d7..02db880d8caaa 100644 --- a/stdlib/SuiteSparse/src/umfpack.jl +++ b/stdlib/SuiteSparse/src/umfpack.jl @@ -387,18 +387,7 @@ function nnz(lu::UmfpackLU) end ### Solve with Factorization -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} = - (lu = translu.parent; ldiv!(B, transpose(lu), copy(B))) -ldiv!(adjlu::Adjoint{T,<:UmfpackLU{T}}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = - (lu = adjlu.parent; ldiv!(B, adjoint(lu), copy(B))) -ldiv!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = - ldiv!(B, lu, copy(B)) -ldiv!(translu::Transpose{Float64,<:UmfpackLU{Float64}}, B::StridedVecOrMat{<:Complex}) = - (lu = translu.parent; ldiv!(B, transpose(lu), copy(B))) -ldiv!(adjlu::Adjoint{Float64,<:UmfpackLU{Float64}}, B::StridedVecOrMat{<:Complex}) = - (lu = adjlu.parent; ldiv!(B, adjoint(lu), copy(B))) +import Base.LinAlg.ldiv! ldiv!(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = _Aq_ldiv_B!(X, lu, B, UMFPACK_A) diff --git a/stdlib/SuiteSparse/test/spqr.jl b/stdlib/SuiteSparse/test/spqr.jl index ff9a93eedf789..59ebeff528bbc 100644 --- a/stdlib/SuiteSparse/test/spqr.jl +++ b/stdlib/SuiteSparse/test/spqr.jl @@ -2,7 +2,7 @@ using SuiteSparse.SPQR using SuiteSparse.CHOLMOD -using Base.LinAlg: mul!, Adjoint, Transpose +using Base.LinAlg: 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 ec747d4a42d1d..e2901cbce3a80 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 = Base.LinAlg.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 = Base.LinAlg.ldiv!(adjoint(lua), z) @test x ≈ float([1:5;]) @test x === z y = similar(x) - SuiteSparse.ldiv!(y, adjoint(lua), complex.(b)) + Base.LinAlg.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 = Base.LinAlg.ldiv!(transpose(lua), complex.(b)) @test x ≈ float([1:5;]) y = similar(x) - SuiteSparse.ldiv!(y, transpose(lua), complex.(b)) + Base.LinAlg.ldiv!(y, transpose(lua), complex.(b)) @test y ≈ x @test transpose(A) * x ≈ b diff --git a/test/linalg/diagonal.jl b/test/linalg/diagonal.jl index ecfd6e6e533bd..bab600a4c9d54 100644 --- a/test/linalg/diagonal.jl +++ b/test/linalg/diagonal.jl @@ -1,8 +1,8 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license using Test, Random -using Base.LinAlg: mul!, ldiv!, rdiv! -import Base.LinAlg: BlasFloat, BlasComplex, SingularException +using Base.LinAlg: mul!, mul1!, mul2!, ldiv!, rdiv!, + BlasFloat, BlasComplex, SingularException using SparseArrays n=12 #Size of matrix problem to test @@ -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 LinAlg.mul1!(copy(D1),D2) == D1*D2 + @test LinAlg.mul2!(D1,copy(D2)) == D1*D2 + @test LinAlg.mul1!(copy(D1),transpose(D2)) == D1*transpose(D2) + @test LinAlg.mul2!(transpose(D1),copy(D2)) == transpose(D1)*D2 + @test LinAlg.mul1!(copy(D1),adjoint(D2)) == D1*adjoint(D2) + @test LinAlg.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/test/linalg/givens.jl b/test/linalg/givens.jl index bb0b04ca9af17..dcf79506b2062 100644 --- a/test/linalg/givens.jl +++ b/test/linalg/givens.jl @@ -1,7 +1,7 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license using Test, Random -using Base.LinAlg: mul! +using Base.LinAlg: mul1!, mul2! # Test givens rotations @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) @@ -16,11 +16,11 @@ using Base.LinAlg: mul! 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) @@ -33,8 +33,8 @@ using Base.LinAlg: mul! @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/test/linalg/lq.jl b/test/linalg/lq.jl index 63e353e7ea407..c201b7eb98728 100644 --- a/test/linalg/lq.jl +++ b/test/linalg/lq.jl @@ -2,7 +2,7 @@ using Test, Random -using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, mul! +using Base.LinAlg: 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::LinAlg.LQPackedQ) = (n = size(Q.factors, 2); mul!(Q, Matrix{eltype(Q)}(I, n, n))) +squareQ(Q::LinAlg.LQPackedQ) = (n = size(Q.factors, 2); mul2!(Q, Matrix{eltype(Q)}(I, n, n))) rectangularQ(Q::LinAlg.LQPackedQ) = convert(Array, Q) @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64) @@ -54,8 +54,6 @@ rectangularQ(Q::LinAlg.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::LinAlg.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::Base.LinAlg.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::Base.LinAlg.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/test/linalg/qr.jl b/test/linalg/qr.jl index 53151902c7988..1b643e9f1ccff 100644 --- a/test/linalg/qr.jl +++ b/test/linalg/qr.jl @@ -2,7 +2,7 @@ using Test, Random -using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, mul! +using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, mul1!, mul2! n = 10 @@ -20,7 +20,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::LinAlg.AbstractQ) = (sq = size(Q.factors, 1); mul!(Q, Matrix{eltype(Q)}(I, sq, sq))) +squareQ(Q::LinAlg.AbstractQ) = (sq = size(Q.factors, 1); mul2!(Q, Matrix{eltype(Q)}(I, sq, sq))) rectangularQ(Q::LinAlg.AbstractQ) = convert(Array, Q) @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) @@ -135,20 +135,20 @@ rectangularQ(Q::LinAlg.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 mu1l!(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 Base.LinAlg.mul!(q,zeros(eltya,n1+1)) - @test_throws DimensionMismatch Base.LinAlg.mul!(adjoint(q), zeros(eltya,n1+1)) + @test_throws DimensionMismatch Base.LinAlg.mul2!(q,zeros(eltya,n1+1)) + @test_throws DimensionMismatch Base.LinAlg.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/test/linalg/special.jl b/test/linalg/special.jl index e63384323d0e9..05092cb864c77 100644 --- a/test/linalg/special.jl +++ b/test/linalg/special.jl @@ -3,7 +3,7 @@ using Test, Random using SparseArrays -using Base.LinAlg: mul! +using Base.LinAlg: mul1! n= 10 #Size of matrix to test srand(1) @@ -100,7 +100,7 @@ end end C = rand(n,n) - for TriType in [Base.LinAlg.UnitLowerTriangular, Base.LinAlg.UnitUpperTriangular, UpperTriangular, LowerTriangular] + for TriType in [LinAlg.UnitLowerTriangular, LinAlg.UnitUpperTriangular, UpperTriangular, LowerTriangular] D = TriType(C) for Spectype in [Diagonal, Bidiagonal, Tridiagonal, Matrix] @test Matrix(D + convert(Spectype,A)) ≈ Matrix(D + A) @@ -112,16 +112,16 @@ end end @testset "Triangular Types and QR" begin - for typ in [UpperTriangular,LowerTriangular,Base.LinAlg.UnitUpperTriangular,Base.LinAlg.UnitLowerTriangular] + for typ in [UpperTriangular,LowerTriangular,LinAlg.UnitUpperTriangular,LinAlg.UnitLowerTriangular] a = rand(n,n) atri = typ(a) 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 @@ -174,8 +174,8 @@ end N = 4 # The tested annotation types testfull = Bool(parse(Int,(get(ENV, "JULIA_TESTFULL", "0")))) - utriannotations = (UpperTriangular, Base.LinAlg.UnitUpperTriangular) - ltriannotations = (LowerTriangular, Base.LinAlg.UnitLowerTriangular) + utriannotations = (UpperTriangular, LinAlg.UnitUpperTriangular) + ltriannotations = (LowerTriangular, LinAlg.UnitLowerTriangular) triannotations = (utriannotations..., ltriannotations...) symannotations = (Symmetric, Hermitian) annotations = testfull ? (triannotations..., symannotations...) : (LowerTriangular, Symmetric) diff --git a/test/linalg/triangular.jl b/test/linalg/triangular.jl index 66bd4b37b6907..8c53fcf9bb60a 100644 --- a/test/linalg/triangular.jl +++ b/test/linalg/triangular.jl @@ -4,7 +4,7 @@ debug = false using Test, Random using Base.LinAlg: BlasFloat, errorbounds, full!, naivesub!, transpose!, UnitUpperTriangular, UnitLowerTriangular, - mul!, rdiv! + mul!, mul1!, mul2!, rdiv! using SparseArrays debug && println("Triangular matrices") @@ -321,7 +321,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 @@ -359,12 +359,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]