From 3e00fb97faa5a0fe810d33854e5fdec0ac397e11 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 6 Jun 2023 22:07:55 +0200 Subject: [PATCH 1/6] Adjust to unwrapping of triangular mul --- src/linalg.jl | 309 ++++++++++++++++++++++++-------------------- src/sparsevector.jl | 2 +- 2 files changed, 173 insertions(+), 138 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index a2f07b22..cbb8597b 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -455,7 +455,7 @@ const TriangularSparse{T} = Union{ LowerTriangularSparse{T}, UpperTriangularSparse{T}} where T ## triangular multipliers -function LinearAlgebra._multrimat!(C::StridedVecOrMat{T}, A::TriangularSparse{T}, B::StridedVecOrMat{T}) where T +function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, tfun::Function, A::SparseMatrixCSCUnion, B::AbstractVecOrMat) C !== B && copyto!(C, B) require_one_based_indexing(A, C) nrowC = size(C, 1) @@ -463,169 +463,203 @@ function LinearAlgebra._multrimat!(C::StridedVecOrMat{T}, A::TriangularSparse{T} if nrowC != ncol throw(DimensionMismatch("A has $(ncol) columns and B has $(nrowC) rows")) end - _lmul!(A, C) -end - -# forward multiplication for UpperTriangular SparseCSC matrices -function _lmul!(U::UpperTriangularPlain, B::StridedVecOrMat) - A = U.data - unit = U isa UnitUpperTriangular - nrowB, ncolB = size(B, 1), size(B, 2) aa = getnzval(A) ja = getrowval(A) ia = getcolptr(A) - joff = 0 - for k = 1:ncolB - for j = 1:nrowB - i1 = ia[j] - i2 = ia[j + 1] - 1 - done = unit - - bj = B[joff + j] - for ii = i1:i2 - jai = ja[ii] - aii = aa[ii] - if jai < j - B[joff + jai] += aii * bj - elseif jai == j - if !unit - B[joff + j] *= aii - done = true + unit = isunitc == 'U' + Z = zero(eltype(C)) + + if uploc == 'U' + if tfun === identity + # forward multiplication for UpperTriangular SparseCSC matrices + for k = 1:ncolB + for j = 1:nrowB + i1 = ia[j] + i2 = ia[j + 1] - 1 + done = unit + + bj = B[joff + j] + for ii = i1:i2 + jai = ja[ii] + aii = aa[ii] + if jai < j + C[joff + jai] += aii * bj + elseif jai == j + if !unit + C[joff + j] = aii * bj + done = true + end + else + break + end + end + if !done + C[joff + j] = Z end - else - break end + joff += nrowB end - if !done - B[joff + j] -= B[joff + j] + else # tfun in (adjoint, transpose) + # backward multiplication with adjoint and transpose of LowerTriangular CSC matrices + for k = 1:ncolB + for j = nrowB:-1:1 + i1 = ia[j] + i2 = ia[j + 1] - 1 + akku = Z + j0 = !unit ? j : j - 1 + + # loop through column j of A - only structural non-zeros + for ii = i1:i2 + jai = ja[ii] + if jai <= j0 + akku += tfun(aa[ii]) * B[joff + jai] + else + break + end + end + if unit + akku += oneunit(eltype(A)) * B[joff + j] + end + C[joff + j] = akku + end + joff += nrowB end end - joff += nrowB - end - B -end - -# backward multiplication for LowerTriangular SparseCSC matrices -function _lmul!(L::LowerTriangularPlain, B::StridedVecOrMat) - A = L.data - unit = L isa UnitLowerTriangular - - nrowB, ncolB = size(B, 1), size(B, 2) - aa = getnzval(A) - ja = getrowval(A) - ia = getcolptr(A) - - joff = 0 - for k = 1:ncolB - for j = nrowB:-1:1 - i1 = ia[j] - i2 = ia[j + 1] - 1 - done = unit - - bj = B[joff + j] - for ii = i2:-1:i1 - jai = ja[ii] - aii = aa[ii] - if jai > j - B[joff + jai] += aii * bj - elseif jai == j - if !unit - B[joff + j] *= aii - done = true + else # uploc == 'L' + if tfun === identity + # backward multiplication for LowerTriangular SparseCSC matrices + for k = 1:ncolB + for j = nrowB:-1:1 + i1 = ia[j] + i2 = ia[j + 1] - 1 + done = unit + + bj = B[joff + j] + for ii = i2:-1:i1 + jai = ja[ii] + aii = aa[ii] + if jai > j + C[joff + jai] += aii * bj + elseif jai == j + if !unit + C[joff + j] = aii * bj + done = true + end + else + break + end + end + if !done + C[joff + j] = Z end - else - break end + joff += nrowB end - if !done - B[joff + j] -= B[joff + j] + else # tfun in (adjoint, transpose) + # forward multiplication for adjoint and transpose of LowerTriangular CSC matrices + for k = 1:ncolB + for j = 1:nrowB + i1 = ia[j] + i2 = ia[j + 1] - 1 + akku = Z + j0 = !unit ? j : j + 1 + + # loop through column j of A - only structural non-zeros + for ii = i2:-1:i1 + jai = ja[ii] + if jai >= j0 + akku += tfun(aa[ii]) * B[joff + jai] + else + break + end + end + if unit + akku += oneunit(eltype(A)) * B[joff + j] + end + C[joff + j] = akku + end + joff += nrowB end end - joff += nrowB end - B + return C end - -# forward multiplication for adjoint and transpose of LowerTriangular CSC matrices -function _lmul!(U::UpperTriangularWrapped, B::StridedVecOrMat) - A = parent(parent(U)) - unit = U isa UnitUpperTriangular - tfun = LinearAlgebra.adj_or_trans(parent(U)) - +function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, _, xA::AdjOrTrans{<:Any,<:SparseMatrixCSCUnion}, B::AbstractVecOrMat) + C !== B && copyto!(C, B) + A = parent(xA) + nrowC = size(C, 1) + ncol = LinearAlgebra.checksquare(A) + if nrowC != ncol + throw(DimensionMismatch("A has $(ncol) columns and B has $(nrowC) rows")) + end nrowB, ncolB = size(B, 1), size(B, 2) aa = getnzval(A) ja = getrowval(A) ia = getcolptr(A) - Z = zero(eltype(A)) - joff = 0 - for k = 1:ncolB - for j = 1:nrowB - i1 = ia[j] - i2 = ia[j + 1] - 1 - akku = Z - j0 = !unit ? j : j + 1 - - # loop through column j of A - only structural non-zeros - for ii = i2:-1:i1 - jai = ja[ii] - if jai >= j0 - aai = tfun(aa[ii]) - akku += B[joff + jai] * aai - else - break + unit = isunitc == 'U' + Z = zero(eltype(C)) + + if uploc == 'U' + for k = 1:ncolB + for j = 1:nrowB + i1 = ia[j] + i2 = ia[j + 1] - 1 + done = unit + + bj = B[joff + j] + for ii = i1:i2 + jai = ja[ii] + aii = conj(aa[ii]) + if jai < j + C[joff + jai] += aii * bj + elseif jai == j + if !unit + C[joff + j] = aii * bj + done = true + end + else + break + end + end + if !done + C[joff + j] = Z end end - if unit - akku += B[joff + j] - end - B[joff + j] = akku + joff += nrowB end - joff += nrowB - end - B -end - -# backward multiplication with adjoint and transpose of LowerTriangular CSC matrices -function _lmul!(L::LowerTriangularWrapped, B::StridedVecOrMat) - A = parent(parent(L)) - unit = L isa UnitLowerTriangular - tfun = LinearAlgebra.adj_or_trans(parent(L)) - - nrowB, ncolB = size(B, 1), size(B, 2) - aa = getnzval(A) - ja = getrowval(A) - ia = getcolptr(A) - Z = zero(eltype(A)) - - joff = 0 - for k = 1:ncolB - for j = nrowB:-1:1 - i1 = ia[j] - i2 = ia[j + 1] - 1 - akku = Z - j0 = !unit ? j : j - 1 - - # loop through column j of A - only structural non-zeros - for ii = i1:i2 - jai = ja[ii] - if jai <= j0 - aai = tfun(aa[ii]) - akku += B[joff + jai] * aai - else - break + else # uploc == 'L' + for k = 1:ncolB + for j = nrowB:-1:1 + i1 = ia[j] + i2 = ia[j + 1] - 1 + done = unit + + bj = B[joff + j] + for ii = i2:-1:i1 + jai = ja[ii] + aii = conj(aa[ii]) + if jai > j + C[joff + jai] += aii * bj + elseif jai == j + if !unit + C[joff + j] = aii * bj + done = true + end + else + break + end + end + if !done + C[joff + j] = Z end end - if unit - akku += B[joff + j] - end - B[joff + j] = akku + joff += nrowB end - joff += nrowB end - B + return C end ## triangular solvers @@ -785,8 +819,9 @@ function LinearAlgebra._ldiv!(C::StridedVector, U::UpperTriangularWrapped, B::St C end -(\)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = ldiv!(L, Array(B)) -#(*)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B)) +(\)(A::Union{UpperTriangular,LowerTriangular}, B::AbstractSparseMatrixCSC) = A \ Array(B) +(\)(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::AbstractSparseMatrixCSC) = A \ Array(B) +# (*)(L::LinearAlgebra.AbstractTriangular, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B)) ## end of triangular diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 6f220a88..08858dce 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1780,7 +1780,7 @@ function (*)(A::_StridedOrTriangularMatrix{Ta}, x::AbstractSparseVector{Tx}) whe length(x) == n || throw(DimensionMismatch()) Ty = promote_op(matprod, eltype(A), eltype(x)) y = Vector{Ty}(undef, m) - mul!(y, A, x, true, false) + mul!(y, A, x) end function LinearAlgebra.generic_matvecmul!(y::AbstractVector, tA, A::StridedMatrix, x::AbstractSparseVector, From 7519e8de6ba80f787a6a81ab3b3df090aa1a6579 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 11 Jun 2023 21:58:51 +0200 Subject: [PATCH 2/6] include triangular solves --- src/linalg.jl | 355 ++++++++++++++++++++++++-------------------- src/sparsevector.jl | 4 +- 2 files changed, 194 insertions(+), 165 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index cbb8597b..31963677 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -424,39 +424,8 @@ function dot(A::AbstractSparseMatrixCSC, B::Union{DenseMatrixUnion,WrapperMatrix end ## triangular sparse handling - -const LowerTriangularPlain{T} = Union{ - LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, - UnitLowerTriangular{T,<:SparseMatrixCSCUnion{T}}} - -const LowerTriangularWrapped{T} = Union{ - LowerTriangular{T,<:Adjoint{T,<:SparseMatrixCSCUnion{T}}}, - UnitLowerTriangular{T,<:Adjoint{T,<:SparseMatrixCSCUnion{T}}}, - LowerTriangular{T,<:Transpose{T,<:SparseMatrixCSCUnion{T}}}, - UnitLowerTriangular{T,<:Transpose{T,<:SparseMatrixCSCUnion{T}}}} where T - -const UpperTriangularPlain{T} = Union{ - UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, - UnitUpperTriangular{T,<:SparseMatrixCSCUnion{T}}} - -const UpperTriangularWrapped{T} = Union{ - UpperTriangular{T,<:Adjoint{T,<:SparseMatrixCSCUnion{T}}}, - UnitUpperTriangular{T,<:Adjoint{T,<:SparseMatrixCSCUnion{T}}}, - UpperTriangular{T,<:Transpose{T,<:SparseMatrixCSCUnion{T}}}, - UnitUpperTriangular{T,<:Transpose{T,<:SparseMatrixCSCUnion{T}}}} where T - -const UpperTriangularSparse{T} = Union{ - UpperTriangularWrapped{T}, UpperTriangularPlain{T}} where T - -const LowerTriangularSparse{T} = Union{ - LowerTriangularWrapped{T}, LowerTriangularPlain{T}} where T - -const TriangularSparse{T} = Union{ - LowerTriangularSparse{T}, UpperTriangularSparse{T}} where T - -## triangular multipliers +## triangular multiplication function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, tfun::Function, A::SparseMatrixCSCUnion, B::AbstractVecOrMat) - C !== B && copyto!(C, B) require_one_based_indexing(A, C) nrowC = size(C, 1) ncol = LinearAlgebra.checksquare(A) @@ -464,6 +433,7 @@ function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, tf throw(DimensionMismatch("A has $(ncol) columns and B has $(nrowC) rows")) end nrowB, ncolB = size(B, 1), size(B, 2) + C !== B && copyto!(C, B) aa = getnzval(A) ja = getrowval(A) ia = getcolptr(A) @@ -587,13 +557,13 @@ function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, tf return C end function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, _, xA::AdjOrTrans{<:Any,<:SparseMatrixCSCUnion}, B::AbstractVecOrMat) - C !== B && copyto!(C, B) A = parent(xA) nrowC = size(C, 1) ncol = LinearAlgebra.checksquare(A) if nrowC != ncol throw(DimensionMismatch("A has $(ncol) columns and B has $(nrowC) rows")) end + C !== B && copyto!(C, B) nrowB, ncolB = size(B, 1), size(B, 2) aa = getnzval(A) ja = getrowval(A) @@ -663,158 +633,217 @@ function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, _, end ## triangular solvers -# forward substitution for LowerTriangular CSC matrices -function LinearAlgebra._ldiv!(C::StridedVector, L::LowerTriangularPlain, B::StridedVector) - A = L.data - unit = L isa UnitLowerTriangular - C !== B && LinearAlgebra._uconvert_copyto!(C, B, oneunit(eltype(L))) - - nrowB = length(B) - aa = getnzval(A) - ja = getrowval(A) - ia = getcolptr(A) - - for j = 1:nrowB - i1 = ia[j] - i2 = ia[j + 1] - one(eltype(ia)) - - # find diagonal element - ii = searchsortedfirst(ja, j, i1, i2, Base.Order.Forward) - jai = ii > i2 ? zero(eltype(ja)) : ja[ii] - - cj = C[j] - # check for zero pivot and divide with pivot - if jai == j - if !unit - cj /= LinearAlgebra._ustrip(aa[ii]) - C[j] = cj - end - ii += 1 - elseif !unit - throw(LinearAlgebra.SingularException(j)) - end - - # update remaining part - for i = ii:i2 - C[ja[i]] -= cj * LinearAlgebra._ustrip(aa[i]) - end +function LinearAlgebra.generic_trimatdiv!(C::StridedVecOrMat, uploc, isunitc, tfun::Function, A::SparseMatrixCSCUnion, B::AbstractVecOrMat) + mA, nA = size(A) + nrowB, ncolB = size(B, 1), size(B, 2) + if nA != nrowB + throw(DimensionMismatch("second dimension of left hand side A, $nA, and first dimension of right hand side B, $nrowB, must be equal")) end - C -end - -# backward substitution for UpperTriangular CSC matrices -function LinearAlgebra._ldiv!(C::StridedVector, U::UpperTriangularPlain, B::StridedVector) - A = U.data - unit = U isa UnitUpperTriangular - C !== B && LinearAlgebra._uconvert_copyto!(C, B, oneunit(eltype(U))) - - nrowB = length(B) + if size(C) != size(B) + throw(DimensionMismatch("size of output, $(size(C)), does not match size of right hand side, $(size(B))")) + end + C !== B && LinearAlgebra._uconvert_copyto!(C, B, oneunit(eltype(A))) aa = getnzval(A) ja = getrowval(A) ia = getcolptr(A) + unit = isunitc == 'U' - for j = nrowB:-1:1 - i1 = ia[j] - i2 = ia[j + 1] - one(eltype(ia)) - - # find diagonal element - ii = searchsortedlast(ja, j, i1, i2, Base.Order.Forward) - jai = ii < i1 ? zero(eltype(ja)) : ja[ii] - - cj = C[j] - # check for zero pivot and divide with pivot - if jai == j - if !unit - cj /= LinearAlgebra._ustrip(aa[ii]) - C[j] = cj - end - ii -= 1 - elseif !unit - throw(LinearAlgebra.SingularException(j)) - end - - # update remaining part - for i = ii:-1:i1 - C[ja[i]] -= cj * LinearAlgebra._ustrip(aa[i]) - end - end - C -end + if uploc == 'L' + if tfun === identity + # forward substitution for LowerTriangular CSC matrices + for k in 1:ncolB + for j = 1:nrowB + i1 = ia[j] + i2 = ia[j + 1] - one(eltype(ia)) -# forward substitution for adjoint and transpose of UpperTriangular CSC matrices -function LinearAlgebra._ldiv!(C::StridedVector, L::LowerTriangularWrapped, B::StridedVector) - A = parent(parent(L)) - unit = L isa UnitLowerTriangular - tfun = LinearAlgebra.adj_or_trans(parent(L)) - C !== B && LinearAlgebra._uconvert_copyto!(C, B, oneunit(eltype(L))) + # find diagonal element + ii = searchsortedfirst(ja, j, i1, i2, Base.Order.Forward) + jai = ii > i2 ? zero(eltype(ja)) : ja[ii] - nrowB = length(B) - aa = getnzval(A) - ja = getrowval(A) - ia = getcolptr(A) + cj = C[j,k] + # check for zero pivot and divide with pivot + if jai == j + if !unit + cj /= LinearAlgebra._ustrip(aa[ii]) + C[j,k] = cj + end + ii += 1 + elseif !unit + throw(LinearAlgebra.SingularException(j)) + end - for j = 1:nrowB - i1 = ia[j] - i2 = ia[j + 1] - 1 - akku = B[j] - done = false - - # loop through column j of A - only structural non-zeros - for ii = i1:i2 - jai = ja[ii] - if jai < j - akku -= C[jai] * tfun(aa[ii]) - elseif jai == j - akku /= unit ? oneunit(eltype(L)) : tfun(aa[ii]) - done = true - break - else - break + # update remaining part + for i = ii:i2 + C[ja[i],k] -= cj * LinearAlgebra._ustrip(aa[i]) + end + end + end + else # tfun in (adjoint, transpose) + # backward substitution for adjoint and transpose of LowerTriangular CSC matrices + for k in 1:ncolB + for j = nrowB:-1:1 + i1 = ia[j] + i2 = ia[j + 1] - 1 + akku = B[j,k] + done = false + + # loop through column j of A - only structural non-zeros + for ii = i2:-1:i1 + jai = ja[ii] + if jai > j + akku -= C[jai,k] * tfun(aa[ii]) + elseif jai == j + akku /= unit ? oneunit(eltype(A)) : tfun(aa[ii]) + done = true + break + else + break + end + end + if !done && !unit + throw(LinearAlgebra.SingularException(j)) + end + C[j,k] = akku + end end end - if !done && !unit - throw(LinearAlgebra.SingularException(j)) + else # uploc == 'U' + if tfun === identity + # backward substitution for UpperTriangular CSC matrices + for k in 1:ncolB + for j = nrowB:-1:1 + i1 = ia[j] + i2 = ia[j + 1] - one(eltype(ia)) + + # find diagonal element + ii = searchsortedlast(ja, j, i1, i2, Base.Order.Forward) + jai = ii < i1 ? zero(eltype(ja)) : ja[ii] + + cj = C[j,k] + # check for zero pivot and divide with pivot + if jai == j + if !unit + cj /= LinearAlgebra._ustrip(aa[ii]) + C[j,k] = cj + end + ii -= 1 + elseif !unit + throw(LinearAlgebra.SingularException(j)) + end + + # update remaining part + for i = ii:-1:i1 + C[ja[i],k] -= cj * LinearAlgebra._ustrip(aa[i]) + end + end + end + else # tfun in (adjoint, transpose) + # forward substitution for adjoint and transpose of UpperTriangular CSC matrices + for k in 1:ncolB + for j = 1:nrowB + i1 = ia[j] + i2 = ia[j + 1] - 1 + akku = B[j,k] + done = false + + # loop through column j of A - only structural non-zeros + for ii = i1:i2 + jai = ja[ii] + if jai < j + akku -= C[jai,k] * tfun(aa[ii]) + elseif jai == j + akku /= unit ? oneunit(eltype(A)) : tfun(aa[ii]) + done = true + break + else + break + end + end + if !done && !unit + throw(LinearAlgebra.SingularException(j)) + end + C[j,k] = akku + end + end end - C[j] = akku end C end +function LinearAlgebra.generic_trimatdiv!(C::StridedVecOrMat, uploc, isunitc, ::Function, xA::AdjOrTrans{<:Any,<:SparseMatrixCSCUnion}, B::AbstractVecOrMat) + A = parent(xA) + mA, nA = size(A) + nrowB, ncolB = size(B, 1), size(B, 2) + if nA != nrowB + throw(DimensionMismatch("second dimension of left hand side A, $nA, and first dimension of right hand side B, $nrowB, must be equal")) + end + if size(C) != size(B) + throw(DimensionMismatch("size of output, $(size(C)), does not match size of right hand side, $(size(B))")) + end + C !== B && LinearAlgebra._uconvert_copyto!(C, B, oneunit(eltype(A))) -# backward substitution for adjoint and transpose of LowerTriangular CSC matrices -function LinearAlgebra._ldiv!(C::StridedVector, U::UpperTriangularWrapped, B::StridedVector) - A = parent(parent(U)) - unit = U isa UnitUpperTriangular - tfun = LinearAlgebra.adj_or_trans(parent(U)) - C !== B && LinearAlgebra._uconvert_copyto!(C, B, oneunit(eltype(U))) - - nrowB = length(B) aa = getnzval(A) ja = getrowval(A) ia = getcolptr(A) + unit = isunitc == 'U' - for j = nrowB:-1:1 - i1 = ia[j] - i2 = ia[j + 1] - 1 - akku = B[j] - done = false - - # loop through column j of A - only structural non-zeros - for ii = i2:-1:i1 - jai = ja[ii] - if jai > j - akku -= C[jai] * tfun(aa[ii]) - elseif jai == j - akku /= unit ? oneunit(eltype(U)) : tfun(aa[ii]) - done = true - break - else - break + if uploc == 'L' + # forward substitution for LowerTriangular CSC matrices + for k in 1:ncolB + for j = 1:nrowB + i1 = ia[j] + i2 = ia[j + 1] - one(eltype(ia)) + + # find diagonal element + ii = searchsortedfirst(ja, j, i1, i2, Base.Order.Forward) + jai = ii > i2 ? zero(eltype(ja)) : ja[ii] + + cj = C[j,k] + # check for zero pivot and divide with pivot + if jai == j + if !unit + cj /= LinearAlgebra._ustrip(conj(aa[ii])) + C[j,k] = cj + end + ii += 1 + elseif !unit + throw(LinearAlgebra.SingularException(j)) + end + + # update remaining part + for i = ii:i2 + C[ja[i],k] -= cj * LinearAlgebra._ustrip(conj(aa[i])) + end end end - if !done && !unit - throw(LinearAlgebra.SingularException(j)) + else # uploc == 'U' + # backward substitution for UpperTriangular CSC matrices + for k in 1:ncolB + for j = nrowB:-1:1 + i1 = ia[j] + i2 = ia[j + 1] - one(eltype(ia)) + + # find diagonal element + ii = searchsortedlast(ja, j, i1, i2, Base.Order.Forward) + jai = ii < i1 ? zero(eltype(ja)) : ja[ii] + + cj = C[j,k] + # check for zero pivot and divide with pivot + if jai == j + if !unit + cj /= LinearAlgebra._ustrip(conj(aa[ii])) + C[j,k] = cj + end + ii -= 1 + elseif !unit + throw(LinearAlgebra.SingularException(j)) + end + + # update remaining part + for i = ii:-1:i1 + C[ja[i],k] -= cj * LinearAlgebra._ustrip(conj(aa[i])) + end + end end - C[j] = akku end C end diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 08858dce..e82fc075 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -4,7 +4,7 @@ import Base: sort!, findall, copy! import LinearAlgebra: promote_to_array_type, promote_to_arrays_ -using LinearAlgebra: adj_or_trans +using LinearAlgebra: wrapperop ### The SparseVector @@ -1981,7 +1981,7 @@ function *(A::AbstractSparseMatrixCSC, x::AbstractSparseVector) end *(xA::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, x::AbstractSparseVector) = - _At_or_Ac_mul_B((a,b) -> adj_or_trans(xA)(a) * b, xA.parent, x, promote_op(matprod, eltype(xA), eltype(x))) + _At_or_Ac_mul_B((a,b) -> wrapperop(xA)(a) * b, xA.parent, x, promote_op(matprod, eltype(xA), eltype(x))) function _At_or_Ac_mul_B(tfun::Function, A::AbstractSparseMatrixCSC{TvA,TiA}, x::AbstractSparseVector{TvX,TiX}, Tv = promote_op(matprod, TvA, TvX)) where {TvA,TiA,TvX,TiX} From 3caefe8f3a1d302de7171071365fb8324af5d47e Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 13 Jun 2023 09:38:14 +0200 Subject: [PATCH 3/6] target eltype in \ --- src/linalg.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 31963677..7b3a9d21 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -848,8 +848,16 @@ function LinearAlgebra.generic_trimatdiv!(C::StridedVecOrMat, uploc, isunitc, :: C end -(\)(A::Union{UpperTriangular,LowerTriangular}, B::AbstractSparseMatrixCSC) = A \ Array(B) -(\)(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::AbstractSparseMatrixCSC) = A \ Array(B) +function (\)(A::Union{UpperTriangular,LowerTriangular}, B::AbstractSparseMatrixCSC) + require_one_based_indexing(B) + TAB = LinearAlgebra._init_eltype(\, eltype(A), eltype(B)) + ldiv!(Matrix{TAB}(undef, size(B)), A, B) +end +function (\)(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::AbstractSparseMatrixCSC) + require_one_based_indexing(B) + TAB = LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B)) + ldiv!(Matrix{TAB}(undef, size(B)), A, B) +end # (*)(L::LinearAlgebra.AbstractTriangular, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B)) ## end of triangular From c602b5ab4517787ae5875ba2bb405d7981b4868c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 13 Jun 2023 15:56:13 +0200 Subject: [PATCH 4/6] improve package load time --- src/linalg.jl | 50 ++++++++++++++++++++++++++------------------------ 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 7b3a9d21..c7f0d1b5 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1,13 +1,11 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license -import LinearAlgebra: checksquare, sym_uplo +using LinearAlgebra: AbstractTriangular, StridedMaybeAdjOrTransMat, checksquare, sym_uplo using Random: rand! # In matrix-vector multiplication, the correct orientation of the vector is assumed. -const DenseMatrixUnion = Union{StridedMatrix, LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, BitMatrix} -const AdjOrTransDenseMatrix = Union{DenseMatrixUnion,Adjoint{<:Any,<:DenseMatrixUnion},Transpose{<:Any,<:DenseMatrixUnion}} +const DenseMatrixUnion = Union{StridedMatrix, BitMatrix} const DenseInputVector = Union{StridedVector, BitVector} -const DenseInputVecOrMat = Union{AdjOrTransDenseMatrix, DenseInputVector} const DenseVecOrMat = Union{DenseMatrixUnion, DenseInputVector} for op ∈ (:+, :-), Wrapper ∈ (:Hermitian, :Symmetric) @@ -29,13 +27,15 @@ for op ∈ (:+, :-) end LinearAlgebra.generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion, B::DenseMatrixUnion, _add::MulAddMul) = - LinearAlgebra._generic_matmatmul!(C, tA, tB, A, B, _add) + spdensemul!(C, tA, tB, A, B, _add) +LinearAlgebra.generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion, B::AbstractTriangular, _add::MulAddMul) = + spdensemul!(C, tA, tB, A, B, _add) LinearAlgebra.generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion, B::DenseInputVector, _add::MulAddMul) = - LinearAlgebra._generic_matmatmul!(C, tA, 'N', A, B, _add) + spdensemul!(C, tA, 'N', A, B, _add) -function LinearAlgebra._generic_matmatmul!(C::StridedVecOrMat, tA, tB, A::SparseMatrixCSCUnion, B::DenseVecOrMat, _add::MulAddMul) +function spdensemul!(C, tA, tB, A, B, _add) if tA == 'N' - _spmul!(C, A, LinearAlgebra.wrap(B, tB), _add.alpha, _add.beta) + _spmatmul!(C, A, LinearAlgebra.wrap(B, tB), _add.alpha, _add.beta) elseif tA == 'T' _At_or_Ac_mul_B!(transpose, C, A, LinearAlgebra.wrap(B, tB), _add.alpha, _add.beta) elseif tA == 'C' @@ -52,7 +52,7 @@ function LinearAlgebra._generic_matmatmul!(C::StridedVecOrMat, tA, tB, A::Sparse return C end -function _spmul!(C::StridedVecOrMat, A::AbstractSparseMatrixCSC, B::DenseInputVecOrMat, α::Number, β::Number) +function _spmatmul!(C, A, B, α, β) size(A, 2) == size(B, 1) || throw(DimensionMismatch()) size(A, 1) == size(C, 1) || throw(DimensionMismatch()) size(B, 2) == size(C, 2) || throw(DimensionMismatch()) @@ -70,12 +70,10 @@ function _spmul!(C::StridedVecOrMat, A::AbstractSparseMatrixCSC, B::DenseInputVe C end -*(A::SparseMatrixCSCUnion{TA}, x::DenseInputVector) where {TA} = - (T = promote_op(matprod, TA, eltype(x)); mul!(similar(x, T, size(A, 1)), A, x, true, false)) -*(A::SparseMatrixCSCUnion{TA}, B::AdjOrTransDenseMatrix) where {TA} = +*(A::SparseMatrixCSCUnion{TA}, B::AbstractTriangular) where {TA} = (T = promote_op(matprod, TA, eltype(B)); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B, true, false)) -function _At_or_Ac_mul_B!(tfun::Function, C::StridedVecOrMat, A::AbstractSparseMatrixCSC, B::DenseInputVecOrMat, α::Number, β::Number) +function _At_or_Ac_mul_B!(tfun::Function, C, A, B, α, β) size(A, 2) == size(C, 1) || throw(DimensionMismatch()) size(A, 1) == size(B, 1) || throw(DimensionMismatch()) size(B, 2) == size(C, 2) || throw(DimensionMismatch()) @@ -94,9 +92,7 @@ function _At_or_Ac_mul_B!(tfun::Function, C::StridedVecOrMat, A::AbstractSparseM C end -*(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, x::DenseInputVector) = - (T = promote_op(matprod, eltype(A), eltype(x)); mul!(similar(x, T, size(A, 1)), A, x, true, false)) -*(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, B::AdjOrTransDenseMatrix) = +*(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, B::AbstractTriangular) = (T = promote_op(matprod, eltype(A), eltype(B)); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B, true, false)) function LinearAlgebra.generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::AbstractSparseMatrixCSC, _add::MulAddMul) @@ -142,10 +138,14 @@ function _spmul!(C::StridedMatrix, X::AdjOrTrans{<:Any,<:DenseMatrixUnion}, A::A end C end -*(X::AdjOrTransDenseMatrix, A::SparseMatrixCSCUnion{TvA}) where {TvA} = +*(X::StridedMaybeAdjOrTransMat, A::SparseMatrixCSCUnion{TvA}) where {TvA} = + (T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) +*(X::Union{BitMatrix,AdjOrTrans{<:Any,BitMatrix}}, A::SparseMatrixCSCUnion{TvA}) where {TvA} = + (T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) +*(X::AbstractTriangular, A::SparseMatrixCSCUnion{TvA}) where {TvA} = (T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) -function _A_mul_Bt_or_Bc!(tfun::Function, C::StridedMatrix, A::AdjOrTransDenseMatrix, B::AbstractSparseMatrixCSC, α::Number, β::Number) +function _A_mul_Bt_or_Bc!(tfun::Function, C::StridedMatrix, A::AbstractMatrix, B::AbstractSparseMatrixCSC, α::Number, β::Number) mA, nA = size(A) nA == size(B, 2) || throw(DimensionMismatch()) mA == size(C, 1) || throw(DimensionMismatch()) @@ -162,10 +162,12 @@ function _A_mul_Bt_or_Bc!(tfun::Function, C::StridedMatrix, A::AdjOrTransDenseMa end C end -*(X::AdjOrTransDenseMatrix, adjA::Adjoint{<:Any,<:AbstractSparseMatrixCSC}) = - (T = promote_op(matprod, eltype(X), eltype(adjA)); mul!(similar(X, T, (size(X, 1), size(adjA, 2))), X, adjA, true, false)) -*(X::AdjOrTransDenseMatrix, tA::Transpose{<:Any,<:AbstractSparseMatrixCSC}) = - (T = promote_op(matprod, eltype(X), eltype(tA)); mul!(similar(X, T, (size(X, 1), size(tA, 2))), X, tA, true, false)) +*(X::StridedMaybeAdjOrTransMat, A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = + (T = promote_op(matprod, eltype(X), eltype(A)); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) +*(X::Union{BitMatrix,AdjOrTrans{<:Any,BitMatrix}}, A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = + (T = promote_op(matprod, eltype(X), eltype(A)); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) +*(X::AbstractTriangular, A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = + (T = promote_op(matprod, eltype(X), eltype(A)); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) # Sparse matrix multiplication as described in [Gustavson, 1978]: # http://dl.acm.org/citation.cfm?id=355796 @@ -428,7 +430,7 @@ end function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, tfun::Function, A::SparseMatrixCSCUnion, B::AbstractVecOrMat) require_one_based_indexing(A, C) nrowC = size(C, 1) - ncol = LinearAlgebra.checksquare(A) + ncol = checksquare(A) if nrowC != ncol throw(DimensionMismatch("A has $(ncol) columns and B has $(nrowC) rows")) end @@ -559,7 +561,7 @@ end function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, _, xA::AdjOrTrans{<:Any,<:SparseMatrixCSCUnion}, B::AbstractVecOrMat) A = parent(xA) nrowC = size(C, 1) - ncol = LinearAlgebra.checksquare(A) + ncol = checksquare(A) if nrowC != ncol throw(DimensionMismatch("A has $(ncol) columns and B has $(nrowC) rows")) end From 7259c916659ba2089a5defa58ba6e811ce6fb30d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 13 Jun 2023 17:08:32 +0200 Subject: [PATCH 5/6] distinguish dense and sparse triangular matrices --- src/linalg.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index c7f0d1b5..85878e6e 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1,10 +1,12 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license -using LinearAlgebra: AbstractTriangular, StridedMaybeAdjOrTransMat, checksquare, sym_uplo +using LinearAlgebra: AbstractTriangular, StridedMaybeAdjOrTransMat, UpperOrLowerTriangular, + checksquare, sym_uplo using Random: rand! # In matrix-vector multiplication, the correct orientation of the vector is assumed. const DenseMatrixUnion = Union{StridedMatrix, BitMatrix} +const DenseTriangular = UpperOrLowerTriangular{<:Any,<:DenseMatrixUnion} const DenseInputVector = Union{StridedVector, BitVector} const DenseVecOrMat = Union{DenseMatrixUnion, DenseInputVector} @@ -70,8 +72,8 @@ function _spmatmul!(C, A, B, α, β) C end -*(A::SparseMatrixCSCUnion{TA}, B::AbstractTriangular) where {TA} = - (T = promote_op(matprod, TA, eltype(B)); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B, true, false)) +*(A::SparseMatrixCSCUnion{TA}, B::DenseTriangular) where {TA} = + (T = promote_op(matprod, TA, eltype(B)); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B)) function _At_or_Ac_mul_B!(tfun::Function, C, A, B, α, β) size(A, 2) == size(C, 1) || throw(DimensionMismatch()) @@ -92,8 +94,8 @@ function _At_or_Ac_mul_B!(tfun::Function, C, A, B, α, β) C end -*(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, B::AbstractTriangular) = - (T = promote_op(matprod, eltype(A), eltype(B)); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B, true, false)) +*(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, B::DenseTriangular) = + (T = promote_op(matprod, eltype(A), eltype(B)); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B)) function LinearAlgebra.generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::AbstractSparseMatrixCSC, _add::MulAddMul) transA = tA == 'N' ? identity : tA == 'T' ? transpose : adjoint @@ -142,7 +144,7 @@ end (T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) *(X::Union{BitMatrix,AdjOrTrans{<:Any,BitMatrix}}, A::SparseMatrixCSCUnion{TvA}) where {TvA} = (T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) -*(X::AbstractTriangular, A::SparseMatrixCSCUnion{TvA}) where {TvA} = +*(X::DenseTriangular, A::SparseMatrixCSCUnion{TvA}) where {TvA} = (T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) function _A_mul_Bt_or_Bc!(tfun::Function, C::StridedMatrix, A::AbstractMatrix, B::AbstractSparseMatrixCSC, α::Number, β::Number) @@ -166,7 +168,7 @@ end (T = promote_op(matprod, eltype(X), eltype(A)); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) *(X::Union{BitMatrix,AdjOrTrans{<:Any,BitMatrix}}, A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = (T = promote_op(matprod, eltype(X), eltype(A)); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) -*(X::AbstractTriangular, A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = +*(X::DenseTriangular, A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = (T = promote_op(matprod, eltype(X), eltype(A)); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false)) # Sparse matrix multiplication as described in [Gustavson, 1978]: @@ -860,7 +862,7 @@ function (\)(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::AbstractSpars TAB = LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B)) ldiv!(Matrix{TAB}(undef, size(B)), A, B) end -# (*)(L::LinearAlgebra.AbstractTriangular, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B)) +# (*)(L::DenseTriangular, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B)) ## end of triangular From 821cddef0f3584840c3394d46d02204b6136317a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 20 Jun 2023 13:18:20 +0200 Subject: [PATCH 6/6] update to LinAlg PR --- src/linalg.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 85878e6e..426572cb 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -637,6 +637,9 @@ function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, _, end ## triangular solvers +_uconvert_copyto!(c, b, oA) = (c .= Ref(oA) .\ b) +_uconvert_copyto!(c::AbstractArray{T}, b::AbstractArray{T}, _) where {T} = copyto!(c, b) + function LinearAlgebra.generic_trimatdiv!(C::StridedVecOrMat, uploc, isunitc, tfun::Function, A::SparseMatrixCSCUnion, B::AbstractVecOrMat) mA, nA = size(A) nrowB, ncolB = size(B, 1), size(B, 2) @@ -646,7 +649,7 @@ function LinearAlgebra.generic_trimatdiv!(C::StridedVecOrMat, uploc, isunitc, tf if size(C) != size(B) throw(DimensionMismatch("size of output, $(size(C)), does not match size of right hand side, $(size(B))")) end - C !== B && LinearAlgebra._uconvert_copyto!(C, B, oneunit(eltype(A))) + C !== B && _uconvert_copyto!(C, B, oneunit(eltype(A))) aa = getnzval(A) ja = getrowval(A) ia = getcolptr(A) @@ -783,7 +786,7 @@ function LinearAlgebra.generic_trimatdiv!(C::StridedVecOrMat, uploc, isunitc, :: if size(C) != size(B) throw(DimensionMismatch("size of output, $(size(C)), does not match size of right hand side, $(size(B))")) end - C !== B && LinearAlgebra._uconvert_copyto!(C, B, oneunit(eltype(A))) + C !== B && _uconvert_copyto!(C, B, oneunit(eltype(A))) aa = getnzval(A) ja = getrowval(A)