From b4d4498a13406453692e916cbffbd0bb0d9e1bc0 Mon Sep 17 00:00:00 2001 From: Jameson Nash Date: Wed, 12 Jul 2023 13:14:59 -0400 Subject: [PATCH 1/6] fully qualify promote_eltypeof call (#411) Missed in https://github.com/JuliaSparse/SparseArrays.jl/pull/407, but causing doctest build to fail when we try to use this. --- src/sparsevector.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 6ee35828..cb7d82b1 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1236,7 +1236,7 @@ function hvcat(rows::Tuple{Vararg{Int}}, X::_SparseConcatGroup...) if anysparse(X...) vcat(_hvcat_rows(rows, X...)...) else - Base.typed_hvcat(promote_eltypeof(X...), rows, X...) + Base.typed_hvcat(Base.promote_eltypeof(X...), rows, X...) end end function _hvcat_rows((row1, rows...)::Tuple{Vararg{Int}}, X::_SparseConcatGroup...) From b4b0e721ada6e8cf5f6391aff4db307be69b0401 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 12 Jul 2023 20:10:38 +0200 Subject: [PATCH 2/6] Intermediate test relaxation (#405) --- test/linalg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/linalg.jl b/test/linalg.jl index cfc9456d..c4de6764 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -138,7 +138,7 @@ end ACa = sparse(trop(AC)) # copied and adjoint @test AT \ B ≈ AC \ B @test ATa \ B ≈ ACa \ B - @test ATa \ sparse(B) == ATa \ B + @test ATa \ sparse(B) ≈ ATa \ B @test Matrix(ATa) \ B ≈ ATa \ B @test ATa * ( ATa \ B ) ≈ B end From 0eb9c0472de0a392afb0bb70aef1cf68ae5bcb2d Mon Sep 17 00:00:00 2001 From: musvaage <112724366+musvaage@users.noreply.github.com> Date: Sat, 15 Jul 2023 11:40:39 -0500 Subject: [PATCH 3/6] fix typos (#414) --- docs/src/index.md | 2 +- src/higherorderfns.jl | 2 +- src/linalg.jl | 2 +- src/solvers/cholmod.jl | 2 +- src/solvers/umfpack.jl | 4 ++-- test/sparsematrix_constructors_indexing.jl | 4 ++-- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 093a3143..5ddae840 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -122,7 +122,7 @@ julia> R = sparsevec(I,V) The inverse of the [`sparse`](@ref) and [`sparsevec`](@ref) functions is [`findnz`](@ref), which retrieves the inputs used to create the sparse array. -[`findall(!iszero, x)`](@ref) returns the cartesian indices of non-zero entries in `x` +[`findall(!iszero, x)`](@ref) returns the Cartesian indices of non-zero entries in `x` (including stored entries equal to zero). ```jldoctest sparse_function diff --git a/src/higherorderfns.jl b/src/higherorderfns.jl index 92597959..5603d3e6 100644 --- a/src/higherorderfns.jl +++ b/src/higherorderfns.jl @@ -111,7 +111,7 @@ const SpBroadcasted2{Style<:SPVM,Axes,F,Args<:Tuple{SparseVecOrMat,SparseVecOrMa # (1) The definitions below provide a common interface to sparse vectors and matrices # sufficient for the purposes of map[!]/broadcast[!]. This interface treats sparse vectors -# as n-by-one sparse matrices which, though technically incorrect, is how broacast[!] views +# as n-by-one sparse matrices which, though technically incorrect, is how broadcast[!] views # sparse vectors in practice. @inline numrows(A::AbstractCompressedVector) = length(A) @inline numrows(A::AbstractSparseMatrixCSC) = size(A, 1) diff --git a/src/linalg.jl b/src/linalg.jl index a2f07b22..2a4c8db4 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1413,7 +1413,7 @@ kron(a::_SparseKronGroup, b::Number) = a * b ## det, inv, cond -inv(A::AbstractSparseMatrixCSC) = error("The inverse of a sparse matrix can often be dense and can cause the computer to run out of memory. If you are sure you have enough memory, please either convert your matrix to a dense matrix, e.g. by calling `Matrix` or if `A` can be factorized, use `\\` on the dense identity matrix, e.g. `A \\ Matrix{eltype(A)}(I, size(A)...)` restrictions of `\\` on sparse lhs applies. Altenatively, `A\\b` is generally preferable to `inv(A)*b`") +inv(A::AbstractSparseMatrixCSC) = error("The inverse of a sparse matrix can often be dense and can cause the computer to run out of memory. If you are sure you have enough memory, please either convert your matrix to a dense matrix, e.g. by calling `Matrix` or if `A` can be factorized, use `\\` on the dense identity matrix, e.g. `A \\ Matrix{eltype(A)}(I, size(A)...)` restrictions of `\\` on sparse lhs applies. Alternatively, `A\\b` is generally preferable to `inv(A)*b`") # TODO diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 075b2fcf..9c832ef7 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -404,7 +404,7 @@ Factor(FC::FactorComponent) = FC.F ################# # Dense wrappers -# The ifelse here may be unecessary. +# The ifelse here may be unnecessary. # nothing different actually occurs in cholmod_l_allocate vs cholmod_allocate AFAICT. # And CHOLMOD has no way of tracking the difference internally (no internal itype field). # This leads me to believe they can be mixed with long and int versions of sparse freely. diff --git a/src/solvers/umfpack.jl b/src/solvers/umfpack.jl index 9f0cd5e8..26d6c8d6 100644 --- a/src/solvers/umfpack.jl +++ b/src/solvers/umfpack.jl @@ -275,7 +275,7 @@ UmfpackWS(F::ATLU, refinement::Bool=has_refinement(F)) = UmfpackWS(F.parent, ref A shallow copy of UmfpackLU to use in multithreaded solve applications. This function duplicates the working space, control, info and lock fields. """ -# Not using simlar helps if the actual needed size has changed as it would need to be resized again +# Not using similar helps if the actual needed size has changed as it would need to be resized again Base.copy(F::UmfpackLU{Tv, Ti}, ws=UmfpackWS(F)) where {Tv, Ti} = UmfpackLU( F.symbolic, @@ -413,7 +413,7 @@ When `check = false`, responsibility for checking the decomposition's validity (via [`issuccess`](@ref)) lies with the user. The permutation `q` can either be a permutation vector or `nothing`. If no permutation vector -is proveded or `q` is `nothing`, UMFPACK's default is used. If the permutation is not zero based, a +is provided or `q` is `nothing`, UMFPACK's default is used. If the permutation is not zero based, a zero based copy is made. See also [`lu`](@ref) diff --git a/test/sparsematrix_constructors_indexing.jl b/test/sparsematrix_constructors_indexing.jl index 562a14c6..1ef19a1d 100644 --- a/test/sparsematrix_constructors_indexing.jl +++ b/test/sparsematrix_constructors_indexing.jl @@ -1077,7 +1077,7 @@ end @test Base.isstored(A, c[1], c[2]) == false end - # `isstored` for adjoint and tranposed matrices: + # `isstored` for adjoint and transposed matrices: for trans in (adjoint, transpose) B = trans(A) stored_indices = [CartesianIndex(j, i) for (j, i) in zip(J, I)] @@ -1390,7 +1390,7 @@ using Base: swaprows!, swapcols! (swaprows!, 2, 3), # Test swapping rows of unequal length (swaprows!, 2, 4), # Test swapping non-adjacent rows (swapcols!, 1, 2), # Test swapping columns where one column is fully sparse - (swapcols!, 2, 3), # Test swapping coulms of unequal length + (swapcols!, 2, 3), # Test swapping columns of unequal length (swapcols!, 2, 4)) # Test swapping non-adjacent columns Scopy = copy(S) Sdense = Array(S) From f8f0f409972b627e05a5e89f58f189f94e3a5666 Mon Sep 17 00:00:00 2001 From: Sobhan Mohammadpour Date: Tue, 18 Jul 2023 14:14:57 -0400 Subject: [PATCH 4/6] bring coverage of fixed SparseMatrixCSC to 100% (#392) Co-authored-by: Daniel Karrasch Co-authored-by: Viral B. Shah --- src/sparsematrix.jl | 2 +- test/fixed.jl | 34 ++++++++++++++++++++++++++++++++-- 2 files changed, 33 insertions(+), 3 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 08d8a21a..2f8bd9ce 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -124,7 +124,7 @@ Get a writable copy of x. See `_unsafe_unfix(x)` """ SparseMatrixCSC(x::FixedSparseCSC) = SparseMatrixCSC(size(x, 1), size(x, 2), copy(parent(getcolptr(x))), - copy(parent(rowval(x))), + copy(parent(rowvals(x))), copy(nonzeros(x))) function sparse_check_Ti(m::Integer, n::Integer, Ti::Type) diff --git a/test/fixed.jl b/test/fixed.jl index 1441929e..59c10346 100644 --- a/test/fixed.jl +++ b/test/fixed.jl @@ -20,6 +20,23 @@ using SparseArrays: AbstractSparseVector, AbstractSparseMatrixCSC, FixedSparseCS @test (r[1] = r[1]; true) end +@testset "SparseMatrixCSC from readonly" begin + + # test that SparseMatrixCSC from readonly does copy + A = sprandn(12, 11, 0.3) + B = SparseMatrixCSC(size(A)..., ReadOnly(getcolptr(A)), ReadOnly(rowvals(A)), nonzeros(A)) + + @test typeof(B) == typeof(A) + @test A == B + + @test getcolptr(A) == getcolptr(B) + @test getcolptr(A) !== getcolptr(B) + @test rowvals(A) == rowvals(B) + @test rowvals(A) !== rowvals(B) + @test nonzeros(A) == nonzeros(B) + @test nonzeros(A) === nonzeros(B) +end + struct_eq(A, B, C...) = struct_eq(A, B) && struct_eq(B, C...) struct_eq(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSC) = getcolptr(A) == getcolptr(B) && rowvals(A) == rowvals(B) @@ -30,6 +47,10 @@ struct_eq(A::AbstractSparseVector, B::AbstractSparseVector) = A = sprandn(10, 10, 0.3) F = FixedSparseCSC(copy(A)) + Ft = FixedSparseCSC{eltype(A),eltype(rowvals(A))}(A) + @test typeof(Ft) == typeof(F) + @test Ft == F + @test struct_eq(F, A) nonzeros(F) .= 0 @test struct_eq(F, A) @@ -69,7 +90,16 @@ struct_eq(A::AbstractSparseVector, B::AbstractSparseVector) = @test typeof(B) == typeof(F) @test struct_eq(B, F) end +@testset "SparseMatrixCSC conversions" begin + A = sprandn(10, 10, 0.3) + F = fixed(copy(A)) + B = SparseMatrixCSC(F) + @test A == B + # fixed(x...) + @test sparse(2I, 3, 3) == sparse(fixed(2I, 3, 3)) + @test SparseArrays._unsafe_unfix(A) == A +end @testset "FixedSparseVector" begin y = sprandn(10, 0.3) x = FixedSparseVector(copy(y)) @@ -88,7 +118,7 @@ end @test f(x, y, z) == 0 t = similar(x) @test typeof(t) == typeof(x) - @test struct_eq(t, x) + @test struct_eq(t, x) end @testset "Issue #190" begin @@ -139,4 +169,4 @@ always_false(x...) = false @test all(iszero, nonzeros(b)) end -end \ No newline at end of file +end From b3872c8b529a4ffd6b1659f71cfcb02a6582fd3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 21 Jul 2023 22:57:23 +0200 Subject: [PATCH 5/6] added warning for iterating while mutating a sparse matrix (#415) --- src/sparsematrix.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 2f8bd9ce..10bd525a 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -294,6 +294,9 @@ column. In conjunction with [`nonzeros`](@ref) and # perform sparse wizardry... end end + +!!! warning + Adding or removing nonzero elements to the matrix may invalidate the `nzrange`, one should not mutate the matrix while iterating. """ nzrange(S::AbstractSparseMatrixCSC, col::Integer) = getcolptr(S)[col]:(getcolptr(S)[col+1]-1) nzrange(S::SparseMatrixCSCView, col::Integer) = nzrange(S.parent, S.indices[2][col]) From cb10c1e3fe74716b6df40876afed4e075b4a6531 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 22 Jul 2023 12:53:47 +0200 Subject: [PATCH 6/6] use unwrapping mechanism for triangular matrices (#396) --- src/linalg.jl | 729 ++++++++++++++++++++++++-------------------- src/sparsevector.jl | 6 +- 2 files changed, 407 insertions(+), 328 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 2a4c8db4..ef56de83 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1,13 +1,13 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license -import LinearAlgebra: 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, LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, BitMatrix} -const AdjOrTransDenseMatrix = Union{DenseMatrixUnion,Adjoint{<:Any,<:DenseMatrixUnion},Transpose{<:Any,<:DenseMatrixUnion}} +const DenseMatrixUnion = Union{StridedMatrix, BitMatrix} +const DenseTriangular = UpperOrLowerTriangular{<:Any,<:DenseMatrixUnion} const DenseInputVector = Union{StridedVector, BitVector} -const DenseInputVecOrMat = Union{AdjOrTransDenseMatrix, DenseInputVector} const DenseVecOrMat = Union{DenseMatrixUnion, DenseInputVector} for op ∈ (:+, :-), Wrapper ∈ (:Hermitian, :Symmetric) @@ -29,13 +29,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 +54,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 +72,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} = - (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::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,10 +94,8 @@ 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) = - (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,10 +140,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::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::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 +164,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::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]: # http://dl.acm.org/citation.cfm?id=355796 @@ -424,369 +428,444 @@ 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 -function LinearAlgebra._multrimat!(C::StridedVecOrMat{T}, A::TriangularSparse{T}, B::StridedVecOrMat{T}) where T - C !== B && copyto!(C, B) +## triangular multiplication +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 - _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) + C !== B && copyto!(C, B) 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) + A = parent(xA) + nrowC = size(C, 1) + ncol = 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) 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 -# 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))) +_uconvert_copyto!(c, b, oA) = (c .= Ref(oA) .\ b) +_uconvert_copyto!(c::AbstractArray{T}, b::AbstractArray{T}, _) where {T} = copyto!(c, b) - 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) - aa = getnzval(A) - ja = getrowval(A) - ia = getcolptr(A) - - 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 + if size(C) != size(B) + throw(DimensionMismatch("size of output, $(size(C)), does not match size of right hand side, $(size(B))")) end - C -end - -# 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))) - - nrowB = length(B) + C !== B && _uconvert_copyto!(C, B, oneunit(eltype(A))) aa = getnzval(A) ja = getrowval(A) ia = getcolptr(A) + unit = isunitc == 'U' + + 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)) + + # 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(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 && _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' + + 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 - 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 + # 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 -(\)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = ldiv!(L, Array(B)) -#(*)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = lmul!(L, 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::DenseTriangular, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B)) ## end of triangular diff --git a/src/sparsevector.jl b/src/sparsevector.jl index cb7d82b1..22b1badb 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 @@ -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, @@ -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}