From 0e76dbb1563644edbc31d361d3f32006d32827fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 11 Oct 2024 14:58:50 +0200 Subject: [PATCH 1/6] added diagonal-sparse multiplication --- src/linalg.jl | 10 ++++++++++ test/linalg.jl | 10 ++++++++++ 2 files changed, 20 insertions(+) diff --git a/src/linalg.jl b/src/linalg.jl index 131a21bc..e9f6595e 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -188,6 +188,16 @@ function _A_mul_Bt_or_Bc!(tfun::Function, C::StridedMatrix, A::AbstractMatrix, B C end +function *(A::Diagonal, b::AbstractSparseVector) + T = promote_eltype(A, b) + nzind_b = nonzeroinds(b) + nzval_b = nonzeros(b) + if isempty(nzind_b) + return zero(T) + end + return sum(A.diag[nzind_b[idx]] * nzval_b[idx] for idx in eachindex(nzind_b)) +end + # Sparse matrix multiplication as described in [Gustavson, 1978]: # http://dl.acm.org/citation.cfm?id=355796 diff --git a/test/linalg.jl b/test/linalg.jl index 45d42d9f..25eab16c 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -673,6 +673,16 @@ end end end +@testset "diagonal - sparse vector mutliplication" begin + for _ in 1:10 + b = spzeros(10) + b[1:3] .= 1:3 + A = Diagonal(randn(10)) + @test norm(A * b - A * Vector(b)) <= 10eps() + @test norm(A * b - Array(A) * b) <= 10eps() + end +end + @testset "sparse matrix * BitArray" begin A = sprand(5,5,0.3) MA = Array(A) From 5c3e46bca5ddaa20e5a9fbed26fae00be72159e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 11 Oct 2024 15:10:56 +0200 Subject: [PATCH 2/6] learning to do maths --- src/linalg.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index e9f6595e..3314a327 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -190,12 +190,13 @@ end function *(A::Diagonal, b::AbstractSparseVector) T = promote_eltype(A, b) + res = similar(b, T) nzind_b = nonzeroinds(b) nzval_b = nonzeros(b) - if isempty(nzind_b) - return zero(T) + for idx in eachindex(nzind_b) + res[nzind_b[idx]] = A.diag[nzind_b[idx]] * nzval_b[idx] end - return sum(A.diag[nzind_b[idx]] * nzval_b[idx] for idx in eachindex(nzind_b)) + return res end # Sparse matrix multiplication as described in [Gustavson, 1978]: From 5dad812cc2472e62c0b5af22da2917fe4c4fdea3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 11 Oct 2024 17:49:22 +0200 Subject: [PATCH 3/6] added Complex test --- test/linalg.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/linalg.jl b/test/linalg.jl index 25eab16c..b2f9b7c4 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -680,6 +680,9 @@ end A = Diagonal(randn(10)) @test norm(A * b - A * Vector(b)) <= 10eps() @test norm(A * b - Array(A) * b) <= 10eps() + Ac = Diagonal(randn(Complex{Float64}, 10)) + @test norm(Ac * b - Ac * Vector(b)) <= 10eps() + @test norm(Ac * b - Array(Ac) * b) <= 10eps() end end From 85e4ca771bd67454b6499ce78c1cd52061c31985 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 15 Oct 2024 21:29:23 +0200 Subject: [PATCH 4/6] added dimension check --- src/linalg.jl | 5 +++++ test/linalg.jl | 2 ++ 2 files changed, 7 insertions(+) diff --git a/src/linalg.jl b/src/linalg.jl index 3314a327..53c7f7a0 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -189,6 +189,11 @@ function _A_mul_Bt_or_Bc!(tfun::Function, C::StridedMatrix, A::AbstractMatrix, B end function *(A::Diagonal, b::AbstractSparseVector) + if size(A, 2) != length(b) + throw( + DimensionMismatch("The dimension of the matrix A $(size(A)) and of the vector b $(length(b))") + ) + end T = promote_eltype(A, b) res = similar(b, T) nzind_b = nonzeroinds(b) diff --git a/test/linalg.jl b/test/linalg.jl index b2f9b7c4..295eba61 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -683,6 +683,8 @@ end Ac = Diagonal(randn(Complex{Float64}, 10)) @test norm(Ac * b - Ac * Vector(b)) <= 10eps() @test norm(Ac * b - Array(Ac) * b) <= 10eps() + @test_throws DimensionMismatch A * [b; 1] + @test_throws DimensionMismatch A * b[1:end-1] end end From 777dff92ac7849e13a8396d4b0ff6dc0b7cd1c06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 24 Oct 2024 18:10:51 +0200 Subject: [PATCH 5/6] Update src/linalg.jl Co-authored-by: Daniel Karrasch --- src/linalg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg.jl b/src/linalg.jl index 53c7f7a0..ae94603f 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -191,7 +191,7 @@ end function *(A::Diagonal, b::AbstractSparseVector) if size(A, 2) != length(b) throw( - DimensionMismatch("The dimension of the matrix A $(size(A)) and of the vector b $(length(b))") + DimensionMismatch(lazy"The dimension of the matrix A $(size(A)) and of the vector b $(length(b))") ) end T = promote_eltype(A, b) From eedb426415c233f385057195282dc6860c30455a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 24 Oct 2024 20:21:30 +0200 Subject: [PATCH 6/6] Update src/linalg.jl Co-authored-by: Daniel Karrasch --- src/linalg.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/linalg.jl b/src/linalg.jl index ae94603f..01e54f71 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -198,8 +198,9 @@ function *(A::Diagonal, b::AbstractSparseVector) res = similar(b, T) nzind_b = nonzeroinds(b) nzval_b = nonzeros(b) + nzval_res = nonzeros(res) for idx in eachindex(nzind_b) - res[nzind_b[idx]] = A.diag[nzind_b[idx]] * nzval_b[idx] + nzval_res[idx] = A.diag[nzind_b[idx]] * nzval_b[idx] end return res end