From dec77b516052725079e4272f4d799bd0096e0d1a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 28 Feb 2024 18:21:28 +0100 Subject: [PATCH 1/3] Adjust matvec and matmatmul! to new internal LinAlg interface --- Project.toml | 2 +- src/SparseArrays.jl | 6 +++--- src/linalg.jl | 38 +++++++++++++++++++------------------- src/sparsevector.jl | 28 ++++++++++++++-------------- 4 files changed, 37 insertions(+), 37 deletions(-) diff --git a/Project.toml b/Project.toml index 4576fc46..051f11da 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "SparseArrays" uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -version = "1.11.0" +version = "1.12.0" [deps] Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index ff554c8e..320f5c0a 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -6,15 +6,15 @@ Support for sparse arrays. Provides `AbstractSparseArray` and subtypes. module SparseArrays using Base: ReshapedArray, promote_op, setindex_shape_check, to_shape, tail, - require_one_based_indexing, promote_eltype + require_one_based_indexing, promote_eltype, @propagate_inbounds, &, | using Base.Order: Forward using LinearAlgebra using LinearAlgebra: AdjOrTrans, AdjointFactorization, TransposeFactorization, matprod, AbstractQ, AdjointQ, HessenbergQ, QRCompactWYQ, QRPackedQ, LQPackedQ, MulAddMul, - UpperOrLowerTriangular + UpperOrLowerTriangular, @stable_muladdmul -import Base: +, -, *, \, /, &, |, xor, ==, zero, @propagate_inbounds +import Base: +, -, *, \, /, ==, zero import LinearAlgebra: mul!, ldiv!, rdiv!, cholesky, adjoint!, diag, eigen, dot, issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, isbanded, cond, diagm, factorize, ishermitian, norm, opnorm, lmul!, rmul!, tril, triu, diff --git a/src/linalg.jl b/src/linalg.jl index ca70ac4c..ac6de59a 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -47,28 +47,28 @@ for op ∈ (:+, :-) end end -generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::DenseMatrixUnion, _add::MulAddMul) = - spdensemul!(C, tA, tB, A, B, _add) -generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::AbstractTriangular, _add::MulAddMul) = - spdensemul!(C, tA, tB, A, B, _add) -generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion2, B::DenseInputVector, _add::MulAddMul) = - spdensemul!(C, tA, 'N', A, B, _add) - -Base.@constprop :aggressive function spdensemul!(C, tA, tB, A, B, _add) +generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::DenseMatrixUnion, alpha::Number, beta::Number) = + spdensemul!(C, tA, tB, A, B, alpha, beta) +generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::AbstractTriangular, alpha::Number, beta::Number) = + spdensemul!(C, tA, tB, A, B, alpha, beta) +generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion2, B::DenseInputVector, alpha::Number, beta::Number) = + spdensemul!(C, tA, 'N', A, B, alpha, beta) + +Base.@constprop :aggressive function spdensemul!(C, tA, tB, A, B, alpha, beta) if tA == 'N' - _spmatmul!(C, A, wrap(B, tB), _add.alpha, _add.beta) + _spmatmul!(C, A, wrap(B, tB), alpha, beta) elseif tA == 'T' - _At_or_Ac_mul_B!(transpose, C, A, wrap(B, tB), _add.alpha, _add.beta) + _At_or_Ac_mul_B!(transpose, C, A, wrap(B, tB), alpha, beta) elseif tA == 'C' - _At_or_Ac_mul_B!(adjoint, C, A, wrap(B, tB), _add.alpha, _add.beta) + _At_or_Ac_mul_B!(adjoint, C, A, wrap(B, tB), alpha, beta) elseif tA in ('S', 's', 'H', 'h') && tB == 'N' rangefun = isuppercase(tA) ? nzrangeup : nzrangelo diagop = tA in ('S', 's') ? identity : real odiagop = tA in ('S', 's') ? transpose : adjoint T = eltype(C) - _mul!(rangefun, diagop, odiagop, C, A, B, T(_add.alpha), T(_add.beta)) + _mul!(rangefun, diagop, odiagop, C, A, B, T(alpha), T(beta)) else - _generic_matmatmul!(C, 'N', 'N', wrap(A, tA), wrap(B, tB), _add) + @stable_muladdmul _generic_matmatmul!(C, 'N', 'N', wrap(A, tA), wrap(B, tB), MulAddMul(alpha, beta)) end return C end @@ -116,14 +116,14 @@ function _At_or_Ac_mul_B!(tfun::Function, C, A, B, α, β) C end -Base.@constprop :aggressive function generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, _add::MulAddMul) +Base.@constprop :aggressive function generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, alpha::Number, beta::Number) transA = tA == 'N' ? identity : tA == 'T' ? transpose : adjoint if tB == 'N' - _spmul!(C, transA(A), B, _add.alpha, _add.beta) + _spmul!(C, transA(A), B, alpha, beta) elseif tB == 'T' - _A_mul_Bt_or_Bc!(transpose, C, transA(A), B, _add.alpha, _add.beta) + _A_mul_Bt_or_Bc!(transpose, C, transA(A), B, alpha, beta) else # tB == 'C' - _A_mul_Bt_or_Bc!(adjoint, C, transA(A), B, _add.alpha, _add.beta) + _A_mul_Bt_or_Bc!(adjoint, C, transA(A), B, alpha, beta) end return C end @@ -319,10 +319,10 @@ function estimate_mulsize(m::Integer, nnzA::Integer, n::Integer, nnzB::Integer, end Base.@constprop :aggressive function generic_matmatmul!(C::SparseMatrixCSCUnion2, tA, tB, A::SparseMatrixCSCUnion2, - B::SparseMatrixCSCUnion2, _add::MulAddMul) + B::SparseMatrixCSCUnion2, alpha::Number, beta::Number) A, tA = tA in ('H', 'h', 'S', 's') ? (wrap(A, tA), 'N') : (A, tA) B, tB = tB in ('H', 'h', 'S', 's') ? (wrap(B, tB), 'N') : (B, tB) - _generic_matmatmul!(C, tA, tB, A, B, _add) + @stable_muladdmul _generic_matmatmul!(C, tA, tB, A, B, MulAddMul(alpha, beta)) end function _generic_matmatmul!(C::SparseMatrixCSCUnion2, tA, tB, A::AbstractVecOrMat, B::AbstractVecOrMat, _add::MulAddMul) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 206a1e06..44f20a57 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1859,28 +1859,28 @@ function (*)(A::_StridedOrTriangularMatrix{Ta}, x::AbstractSparseVector{Tx}) whe end Base.@constprop :aggressive function generic_matvecmul!(y::AbstractVector, tA, A::StridedMatrix, x::AbstractSparseVector, - _add::MulAddMul = MulAddMul()) + alpha::Number, beta::Number) if tA == 'N' - _spmul!(y, A, x, _add.alpha, _add.beta) + _spmul!(y, A, x, alpha, beta) elseif tA == 'T' - _At_or_Ac_mul_B!(transpose, y, A, x, _add.alpha, _add.beta) + _At_or_Ac_mul_B!(transpose, y, A, x, alpha, beta) elseif tA == 'C' - _At_or_Ac_mul_B!(adjoint, y, A, x, _add.alpha, _add.beta) + _At_or_Ac_mul_B!(adjoint, y, A, x, alpha, beta) else - _spmul!(y, wrap(A, tA), x, _add.alpha, _add.beta) + _spmul!(y, wrap(A, tA), x, alpha, beta) end return y end function generic_matvecmul!(y::AbstractVector, tA, A::UpperOrLowerTriangular, x::AbstractSparseVector, - _add::MulAddMul = MulAddMul()) + alpha::Number, beta::Number) @assert tA == 'N' Adata = parent(A) if Adata isa Transpose - _At_or_Ac_mul_B!(transpose, y, _fliptri(A), x, _add.alpha, _add.beta) + _At_or_Ac_mul_B!(transpose, y, _fliptri(A), x, alpha, beta) elseif Adata isa Adjoint - _At_or_Ac_mul_B!(adjoint, y, _fliptri(A), x, _add.alpha, _add.beta) + _At_or_Ac_mul_B!(adjoint, y, _fliptri(A), x, alpha, beta) else # Adata is plain - _spmul!(y, A, x, _add.alpha, _add.beta) + _spmul!(y, A, x, alpha, beta) end return y end @@ -1990,15 +1990,15 @@ end # * and mul! Base.@constprop :aggressive function generic_matvecmul!(y::AbstractVector, tA, A::AbstractSparseMatrixCSC, x::AbstractSparseVector, - _add::MulAddMul = MulAddMul()) + alpha::Number, beta::Number) if tA == 'N' - _spmul!(y, A, x, _add.alpha, _add.beta) + _spmul!(y, A, x, alpha, beta) elseif tA == 'T' - _At_or_Ac_mul_B!((a,b) -> transpose(a) * b, y, A, x, _add.alpha, _add.beta) + _At_or_Ac_mul_B!((a,b) -> transpose(a) * b, y, A, x, alpha, beta) elseif tA == 'C' - _At_or_Ac_mul_B!((a,b) -> adjoint(a) * b, y, A, x, _add.alpha, _add.beta) + _At_or_Ac_mul_B!((a,b) -> adjoint(a) * b, y, A, x, alpha, beta) else - LinearAlgebra._generic_matvecmul!(y, 'N', wrap(A, tA), x, _add) + @stable_muladdmul LinearAlgebra._generic_matvecmul!(y, 'N', wrap(A, tA), x, MulAddMul(alpha, beta)) end return y end From 26c60d96b00a33bd8f9c16a447e0a783d939c5b3 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 2 Mar 2024 11:46:38 +0100 Subject: [PATCH 2/3] make changes non-breaking --- src/linalg.jl | 18 +++++++++++++++++- src/sparsevector.jl | 11 +++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/src/linalg.jl b/src/linalg.jl index ac6de59a..1ce09b63 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -53,6 +53,13 @@ generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::Abstra spdensemul!(C, tA, tB, A, B, alpha, beta) generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion2, B::DenseInputVector, alpha::Number, beta::Number) = spdensemul!(C, tA, 'N', A, B, alpha, beta) +# legacy methods: TODO: remove +generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::DenseMatrixUnion, _add::MulAddMul) = + spdensemul!(C, tA, tB, A, B, _add.alpha, _add.beta) +generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::AbstractTriangular, _add::MulAddMul) = + spdensemul!(C, tA, tB, A, B, _add.alpha, _add.beta) +generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion2, B::DenseInputVector, _add::MulAddMul) = + spdensemul!(C, tA, 'N', A, B, _add.alpha, _add.beta) Base.@constprop :aggressive function spdensemul!(C, tA, tB, A, B, alpha, beta) if tA == 'N' @@ -68,7 +75,7 @@ Base.@constprop :aggressive function spdensemul!(C, tA, tB, A, B, alpha, beta) T = eltype(C) _mul!(rangefun, diagop, odiagop, C, A, B, T(alpha), T(beta)) else - @stable_muladdmul _generic_matmatmul!(C, 'N', 'N', wrap(A, tA), wrap(B, tB), MulAddMul(alpha, beta)) + @stable_muladdmul LinearAlgebra._generic_matmatmul!(C, 'N', 'N', wrap(A, tA), wrap(B, tB), MulAddMul(alpha, beta)) end return C end @@ -116,6 +123,9 @@ function _At_or_Ac_mul_B!(tfun::Function, C, A, B, α, β) C end +# TODO:remove +generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, _add::MulAddMul) = + generic_matmatmul!(C, tA, tB, A, B, _add.alpha, _add.beta) Base.@constprop :aggressive function generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, alpha::Number, beta::Number) transA = tA == 'N' ? identity : tA == 'T' ? transpose : adjoint if tB == 'N' @@ -318,6 +328,12 @@ function estimate_mulsize(m::Integer, nnzA::Integer, n::Integer, nnzB::Integer, p >= 1 ? m*k : p > 0 ? Int(ceil(-expm1(log1p(-p) * n)*m*k)) : 0 # (1-(1-p)^n)*m*k end +# TODO: remove this one method +Base.@constprop :aggressive function generic_matmatmul!(C::SparseMatrixCSCUnion2, tA, tB, A::SparseMatrixCSCUnion2, B::SparseMatrixCSCUnion2, _add::MulAddMul) + A, tA = tA in ('H', 'h', 'S', 's') ? (wrap(A, tA), 'N') : (A, tA) + B, tB = tB in ('H', 'h', 'S', 's') ? (wrap(B, tB), 'N') : (B, tB) + _generic_matmatmul!(C, tA, tB, A, B, _add) +end Base.@constprop :aggressive function generic_matmatmul!(C::SparseMatrixCSCUnion2, tA, tB, A::SparseMatrixCSCUnion2, B::SparseMatrixCSCUnion2, alpha::Number, beta::Number) A, tA = tA in ('H', 'h', 'S', 's') ? (wrap(A, tA), 'N') : (A, tA) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 44f20a57..b77f901d 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1858,6 +1858,10 @@ function (*)(A::_StridedOrTriangularMatrix{Ta}, x::AbstractSparseVector{Tx}) whe mul!(y, A, x) end +# TODO: remove +Base.@constprop :aggressive generic_matvecmul!(y::AbstractVector, tA, A::StridedMatrix, x::AbstractSparseVector, + _add::MulAddMul = MulAddMul()) = + generic_matvecmul!(y, tA, A, x, _add.alpha, _add.beta) Base.@constprop :aggressive function generic_matvecmul!(y::AbstractVector, tA, A::StridedMatrix, x::AbstractSparseVector, alpha::Number, beta::Number) if tA == 'N' @@ -1871,6 +1875,9 @@ Base.@constprop :aggressive function generic_matvecmul!(y::AbstractVector, tA, A end return y end +# TODO: remove +generic_matvecmul!(y::AbstractVector, tA, A::UpperOrLowerTriangular, x::AbstractSparseVector, _add::MulAddMul = MulAddMul()) = + generic_matvecmul!(y, tA, A, x, _add.alpha, _add.beta) function generic_matvecmul!(y::AbstractVector, tA, A::UpperOrLowerTriangular, x::AbstractSparseVector, alpha::Number, beta::Number) @assert tA == 'N' @@ -1989,6 +1996,10 @@ function densemv(A::AbstractSparseMatrixCSC, x::AbstractSparseVector; trans::Abs end # * and mul! +# TODO: remove +Base.@constprop :aggressive generic_matvecmul!(y::AbstractVector, tA, A::AbstractSparseMatrixCSC, x::AbstractSparseVector, + _add::MulAddMul = MulAddMul()) = + generic_matvecmul!(y, tA, A, x, _add.alpha, _add.beta) Base.@constprop :aggressive function generic_matvecmul!(y::AbstractVector, tA, A::AbstractSparseMatrixCSC, x::AbstractSparseVector, alpha::Number, beta::Number) if tA == 'N' From 48f3784d985f5774132be472809e46c5183e5a45 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 8 May 2024 11:06:19 +0200 Subject: [PATCH 3/3] further adjustments --- src/linalg.jl | 36 +++++++++++------------------------- 1 file changed, 11 insertions(+), 25 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 1ce09b63..131a21bc 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -53,25 +53,19 @@ generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::Abstra spdensemul!(C, tA, tB, A, B, alpha, beta) generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion2, B::DenseInputVector, alpha::Number, beta::Number) = spdensemul!(C, tA, 'N', A, B, alpha, beta) -# legacy methods: TODO: remove -generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::DenseMatrixUnion, _add::MulAddMul) = - spdensemul!(C, tA, tB, A, B, _add.alpha, _add.beta) -generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::AbstractTriangular, _add::MulAddMul) = - spdensemul!(C, tA, tB, A, B, _add.alpha, _add.beta) -generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion2, B::DenseInputVector, _add::MulAddMul) = - spdensemul!(C, tA, 'N', A, B, _add.alpha, _add.beta) Base.@constprop :aggressive function spdensemul!(C, tA, tB, A, B, alpha, beta) - if tA == 'N' + tA_uc, tB_uc = uppercase(tA), uppercase(tB) + if tA_uc == 'N' _spmatmul!(C, A, wrap(B, tB), alpha, beta) - elseif tA == 'T' + elseif tA_uc == 'T' _At_or_Ac_mul_B!(transpose, C, A, wrap(B, tB), alpha, beta) - elseif tA == 'C' + elseif tA_uc == 'C' _At_or_Ac_mul_B!(adjoint, C, A, wrap(B, tB), alpha, beta) - elseif tA in ('S', 's', 'H', 'h') && tB == 'N' + elseif tA_uc in ('S', 'H') && tB_uc == 'N' rangefun = isuppercase(tA) ? nzrangeup : nzrangelo - diagop = tA in ('S', 's') ? identity : real - odiagop = tA in ('S', 's') ? transpose : adjoint + diagop = tA_uc == 'S' ? identity : real + odiagop = tA_uc == 'S' ? transpose : adjoint T = eltype(C) _mul!(rangefun, diagop, odiagop, C, A, B, T(alpha), T(beta)) else @@ -123,9 +117,6 @@ function _At_or_Ac_mul_B!(tfun::Function, C, A, B, α, β) C end -# TODO:remove -generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, _add::MulAddMul) = - generic_matmatmul!(C, tA, tB, A, B, _add.alpha, _add.beta) Base.@constprop :aggressive function generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, alpha::Number, beta::Number) transA = tA == 'N' ? identity : tA == 'T' ? transpose : adjoint if tB == 'N' @@ -328,17 +319,12 @@ function estimate_mulsize(m::Integer, nnzA::Integer, n::Integer, nnzB::Integer, p >= 1 ? m*k : p > 0 ? Int(ceil(-expm1(log1p(-p) * n)*m*k)) : 0 # (1-(1-p)^n)*m*k end -# TODO: remove this one method -Base.@constprop :aggressive function generic_matmatmul!(C::SparseMatrixCSCUnion2, tA, tB, A::SparseMatrixCSCUnion2, B::SparseMatrixCSCUnion2, _add::MulAddMul) - A, tA = tA in ('H', 'h', 'S', 's') ? (wrap(A, tA), 'N') : (A, tA) - B, tB = tB in ('H', 'h', 'S', 's') ? (wrap(B, tB), 'N') : (B, tB) - _generic_matmatmul!(C, tA, tB, A, B, _add) -end Base.@constprop :aggressive function generic_matmatmul!(C::SparseMatrixCSCUnion2, tA, tB, A::SparseMatrixCSCUnion2, B::SparseMatrixCSCUnion2, alpha::Number, beta::Number) - A, tA = tA in ('H', 'h', 'S', 's') ? (wrap(A, tA), 'N') : (A, tA) - B, tB = tB in ('H', 'h', 'S', 's') ? (wrap(B, tB), 'N') : (B, tB) - @stable_muladdmul _generic_matmatmul!(C, tA, tB, A, B, MulAddMul(alpha, beta)) + tA_uc, tB_uc = uppercase(tA), uppercase(tB) + Anew, ta = tA_uc in ('S', 'H') ? (wrap(A, tA), oftype(tA, 'N')) : (A, tA) + Bnew, tb = tB_uc in ('S', 'H') ? (wrap(B, tB), oftype(tB, 'N')) : (B, tB) + @stable_muladdmul _generic_matmatmul!(C, ta, tb, Anew, Bnew, MulAddMul(alpha, beta)) end function _generic_matmatmul!(C::SparseMatrixCSCUnion2, tA, tB, A::AbstractVecOrMat, B::AbstractVecOrMat, _add::MulAddMul)