Write tests and fix bug.

Author: RomeoV

src/ThreadedDenseSparseMul.jl

@@ -47,7 +47,8 @@ Fast outer product when computing `C .= β*C + α * a*b'`, but way faster than B
 Also see `fastdensesparse_outer_threaded!` for a multi-threaded version using `Polyester.jl`.
 """
 function fastdensesparse_outer!(C::MatOrView{T}, a::VecOrView{T}, b::SparseVector{T}, α::Number, β::Number) where T
-    C[:, nonzeroinds(b)] .+=  a * nonzeros(b)'
+    C[:, nonzeroinds(b)] .*=  β
+    C[:, nonzeroinds(b)] .+=  a * (α.*nonzeros(b))'
     return C
 end
 

test/runtests.jl

@@ -3,46 +3,59 @@ using Test
 import SparseArrays: mul!, sprand
 import Profile
 
-∈ₛ(a,b) = occursin(a,b)
-# We use a hack to find out if Polyester is actually used, namely
-# profiling the call and checking the used modules for Polyester.
-# This only works if the problem size is large enough (and then is still a statistical property).
-# I wasn't able to find a clear way, but probably there is one.
-function check_module_used(mod_string::String, prof_data::Vector)
-  cache = Dict{Symbol, String}();
-  sframes = Profile.getdict(prof_data) |> values |> Iterators.flatten
-  used_paths = [Profile.short_path(frame.file, cache) for frame in sframes]
-  any(path->mod_string ∈ₛ path, used_paths)
-end
 
 @testset "Compare with equivalent dense mul" begin
-  @testset for trial in 1:10
-    lhs = rand(500, 1000);
-    rhs = sprand(1000, 10_000, 0.1);
+  @testset "fastdensesparse" begin
+    @testset for trial in 1:10
+      lhs = rand(500, 1000);
+      rhs = sprand(1000, 10_000, 0.1);
+
+      baseline = lhs * Matrix(rhs);
+
+      buf = similar(baseline)
+
+      fastdensesparse!(buf, lhs, rhs, 1, 0)
+      @test buf ≈ baseline
+      fastdensesparse!(buf, lhs, rhs, -1, 0)
+      fastdensesparse_threaded!(buf, lhs, rhs, 1, 0)
+      @test buf ≈ baseline
 
-    baseline = lhs * Matrix(rhs);
-    @test lhs * rhs ≈ baseline
 
-    buf = similar(baseline)
+      # test @view interface and \beta \neq 0
+      inds = collect(3:5:100)
+      baseline[inds, :] .+= 2.5 * @view(lhs[inds, :]) * Matrix(rhs);
 
-    @testset "Check actually using Polyester" begin
-        # make sure I've overwritten the regular mul correctly and actually use Polyester...
-        Profile.clear(); Profile.@profile (buf .= lhs * rhs)
-        @test check_module_used("@Polyester", Profile.fetch())
+      fastdensesparse!(@view(buf[inds, :]), @view(lhs[inds, :]), rhs, 2.5, 1)
+      @test buf ≈ baseline rtol=sqrt(eps(eltype(baseline)))
+      fastdensesparse!(@view(buf[inds, :]), @view(lhs[inds, :]), rhs, -2.5, 1)
+      fastdensesparse_threaded!(@view(buf[inds, :]), @view(lhs[inds, :]), rhs, 2.5, 1)
+      @test buf ≈ baseline rtol=sqrt(eps(eltype(baseline)))
     end
+  end
+
+  @testset "fastdensesparse_outer" begin
+    @testset for trial in 1:10
+      lhs = rand(500, 1000);
+      rhs = sprand(1000, 10_000, 0.1);
+      k = rand(1:size(rhs, 1))
 
-    buf .= lhs * rhs
-    @test buf ≈ baseline
+      baseline = lhs[:, k:k] * Matrix(rhs)[k:k, :];
 
-    buf .= 0.
-    buf .+= lhs * rhs
-    @test buf ≈ baseline
+      buf = similar(baseline)
 
-    buf .= 0.
-    buf .+= 2.0.*lhs * rhs
-    @test buf/2 ≈ baseline
+      fastdensesparse_outer!(buf, @view(lhs[:, k]), rhs[k, :], 1, 0)
+      @test buf ≈ baseline
+      fastdensesparse_outer!(buf, @view(lhs[:, k]), rhs[k, :], -1, 0)
+      fastdensesparse_outer_threaded!(buf, @view(lhs[:, k]), rhs[k, :], 1, 0)
+      @test buf ≈ baseline
 
-    mul!(buf, lhs, rhs, 1, 0)
-    @test buf ≈ baseline
+      baseline .+= 2.5 * lhs[:, (k+1):(k+1)] * Matrix(rhs)[(k+1):(k+1), :];
+
+      fastdensesparse_outer!(buf, lhs[:, k+1], rhs[k+1, :], 2.5, 1)
+      @test buf ≈ baseline
+      fastdensesparse_outer!(buf, lhs[:, k+1], rhs[k+1, :], -2.5, 1)
+      fastdensesparse_outer_threaded!(buf, lhs[:, k+1], rhs[k+1, :], 2.5, 1)
+      @test buf ≈ baseline
+    end
   end
 end