diff --git a/Project.toml b/Project.toml
index 2abfd48..d93698d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SparseArraysBase"
 uuid = "0d5efcca-f356-4864-8770-e1ed8d78f208"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.2.7"
+version = "0.2.8"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
diff --git a/src/abstractsparsearrayinterface.jl b/src/abstractsparsearrayinterface.jl
index 65937f3..6d83385 100644
--- a/src/abstractsparsearrayinterface.jl
+++ b/src/abstractsparsearrayinterface.jl
@@ -357,11 +357,15 @@ function sparse_mul!(
   β::Number=false;
   (mul!!)=(default_mul!!),
 )
+  # TODO: Change to: `a_dest .*= β`
+  # once https://github.com/ITensor/SparseArraysBase.jl/issues/19 is fixed.
+  storedvalues(a_dest) .*= β
+  β′ = one(Bool)
   for I1 in eachstoredindex(a1)
     for I2 in eachstoredindex(a2)
       I_dest = mul_indices(I1, I2)
       if !isnothing(I_dest)
-        a_dest[I_dest] = mul!!(a_dest[I_dest], a1[I1], a2[I2], α, β)
+        a_dest[I_dest] = mul!!(a_dest[I_dest], a1[I1], a2[I2], α, β′)
       end
     end
   end
diff --git a/test/Project.toml b/test/Project.toml
index 4e5feb3..575c9b2 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -5,7 +5,10 @@ ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 DerivableInterfaces = "6c5e35bf-e59e-4898-b73c-732dcc4ba65f"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/basics/test_linalg.jl b/test/basics/test_linalg.jl
new file mode 100644
index 0000000..f48585c
--- /dev/null
+++ b/test/basics/test_linalg.jl
@@ -0,0 +1,38 @@
+using SparseArraysBase: SparseArrayDOK
+using LinearAlgebra: mul!
+using Random: Random
+using StableRNGs: StableRNG
+
+const rng = StableRNG(123)
+
+# TODO: add this to main package
+function sprand(rng::Random.AbstractRNG, ::Type{T}, sz::Base.Dims; p::Real=0.5) where {T}
+  A = SparseArrayDOK{T}(undef, sz)
+  for I in eachindex(A)
+    if rand(rng) < p
+      A[I] = rand(rng, T)
+    end
+  end
+  return A
+end
+
+@testset "mul!" begin
+  T = Float64
+  szA = (2, 2)
+  szB = (2, 2)
+  szC = (szA[1], szB[2])
+
+  for p in 0.0:0.25:1
+    C = sprand(rng, T, szC; p)
+    A = sprand(rng, T, szA; p)
+    B = sprand(rng, T, szB; p)
+
+    check1 = mul!(Array(C), Array(A), Array(B))
+    @test mul!(copy(C), A, B) ≈ check1
+
+    α = rand(rng, T)
+    β = rand(rng, T)
+    check2 = mul!(Array(C), Array(A), Array(B), α, β)
+    @test mul!(copy(C), A, B, α, β) ≈ check2
+  end
+end