From 990d37758ffb92520ba501c1db8baf173c31c5dc Mon Sep 17 00:00:00 2001
From: maltezfaria <maltezfaria@gmail.com>
Date: Thu, 10 Apr 2025 09:39:47 +0200
Subject: [PATCH] fix bug when adding sparse to `HMatrix` - previously errored
 when the `eltype` was not a scalar

---
 src/hmatrix.jl       |  2 +-
 test/hmatrix_test.jl | 65 ++++++++++++++++++++++++++++++++++++++------
 2 files changed, 57 insertions(+), 10 deletions(-)

diff --git a/src/hmatrix.jl b/src/hmatrix.jl
index de5a0e0..4d39ef3 100644
--- a/src/hmatrix.jl
+++ b/src/hmatrix.jl
@@ -567,7 +567,7 @@ function _axpy!(a, X::AbstractSparseArray, Y::HMatrix)
                         a = zeros(T, m)
                         b = zeros(T, n)
                         a[i-irange.start+1] = vals[idx]
-                        b[j-jrange.start+1] = 1
+                        b[j-jrange.start+1] = one(T)
                         R.A = hcat(R.A, a)
                         R.B = hcat(R.B, b)
                     else
diff --git a/test/hmatrix_test.jl b/test/hmatrix_test.jl
index a8e0cf3..ccdddfe 100644
--- a/test/hmatrix_test.jl
+++ b/test/hmatrix_test.jl
@@ -37,7 +37,7 @@ include(joinpath(HMatrices.PROJECT_ROOT, "test", "testutils.jl"))
     end
 end
 
-@testset "Sparse arrays" begin
+@testset "Band sparse arrays" begin
     m = 1000
     n = 1000
 
@@ -49,12 +49,59 @@ end
     adm = StrongAdmissibilityStd(; eta = 3)
     rtol = 1e-5
     comp = PartialACA(; rtol = rtol)
-    K = laplace_matrix(X, Y)
-    H = assemble_hmatrix(K, Xclt, Yclt; adm, comp, threads = false, distributed = false)
-    H_full = Matrix(H)
-    T = eltype(H)
-    m, n = size(H)
-    S = spdiagm(0 => rand(T, n))
-    Hnew = axpy!(true, S, deepcopy(H))
-    @test Matrix(Hnew) == (H_full + Matrix(S))
+    @testset "Scalar problem (laplace)" begin
+        K = laplace_matrix(X, Y)
+        H = assemble_hmatrix(K, Xclt, Yclt; adm, comp, threads = false, distributed = false)
+        H_full = Matrix(H)
+        T = eltype(H)
+        m, n = size(H)
+        S = spdiagm(0 => rand(T, n))
+        Hnew = axpy!(true, S, deepcopy(H))
+        @test Matrix(Hnew) == (H_full + Matrix(S))
+    end
+    @testset "Vector problem (elasticity)" begin
+        K = elastostatic_matrix(X, Y, 1, 1)
+        H = assemble_hmatrix(K, Xclt, Yclt; adm, comp, threads = false, distributed = false)
+        H_full = Matrix(H)
+        T = eltype(H)
+        m, n = size(H)
+        S = spdiagm(0 => rand(T, n))
+        Hnew = axpy!(true, S, deepcopy(H))
+        @test Matrix(Hnew) == (H_full + Matrix(S))
+    end
+end
+
+@testset "Precision of adding any sparse matrix to an H-matrix with regard to the matrix - vector multiplication" begin
+    m = 1000
+    n = 1000
+
+    X = rand(SVector{3,Float64}, m)
+    Y = X
+    splitter = CardinalitySplitter(; nmax = 40)
+    Xclt = ClusterTree(X, splitter)
+    Yclt = ClusterTree(Y, splitter)
+    adm = StrongAdmissibilityStd(; eta = 3)
+    rtol = 1e-5
+    comp = PartialACA(; rtol = rtol)
+    U = rand(SVector{3,Float64}, m)
+    @testset "Scalar problem (laplace)" begin
+        K = laplace_matrix(X, Y)
+        H = assemble_hmatrix(K, Xclt, Yclt; adm, comp, threads = false, distributed = false)
+        H_full = Matrix(H)
+        T = eltype(H)
+        m, n = size(H)
+        S = sprand(T, m, n, 0.01)
+        Hnew = axpy!(true, S, deepcopy(H))
+        @test norm(Hnew * U - (H_full + Matrix(S)) * U) < 1e-5 * norm(Hnew * U)
+    end
+    @testset "Vector problem (elasticity)" begin
+        K = elastostatic_matrix(X, Y, 1, 1)
+        H = assemble_hmatrix(K, Xclt, Yclt; adm, comp, threads = false, distributed = false)
+        H_full = Matrix(H)
+        T = eltype(H)
+        m, n = size(H)
+        S = sprand(T, m, n, 0.1)
+        Hnew = axpy!(true, S, deepcopy(H))
+        @test norm(Hnew * U - (H_full + Matrix(S)) * U) < 1e-5 * norm(Hnew * U)
+    end
 end