From de7a2e2154dc9e30bbc79435aa8d954a1050b3b2 Mon Sep 17 00:00:00 2001
From: Joshua <joshua.tetzner@uni-rostock.de>
Date: Tue, 9 Sep 2025 10:16:47 +0200
Subject: [PATCH 1/2] BEAST extension

I wrote an extensions allowing to use HMatrices for the boundary element
toolkit BEAST.jl.
The extension adds BEAST.jl as a weak dependency which means that BEAST.jl is not
required as a packege but an optional add on. Only if BEAST.jl is loaded before
HMatrices the extension will be available.
---
 Project.toml                 |  7 ++++
 ext/HBEAST/HBEAST.jl         |  9 ++++
 ext/HBEAST/abstractmatrix.jl | 79 ++++++++++++++++++++++++++++++++++++
 src/HMatrices.jl             |  4 ++
 test/Project.toml            |  2 +
 test/extBEAST_test.jl        | 43 ++++++++++++++++++++
 test/runtests.jl             |  3 ++
 7 files changed, 147 insertions(+)
 create mode 100644 ext/HBEAST/HBEAST.jl
 create mode 100644 ext/HBEAST/abstractmatrix.jl
 create mode 100644 test/extBEAST_test.jl

diff --git a/Project.toml b/Project.toml
index 0533e29..59107bb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,7 +11,14 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
+[weakdeps]
+BEAST = "bb4162c7-ba94-5a20-af32-d8ec4428bdd1"
+
+[extensions]
+HBEAST = ["BEAST"]
+
 [compat]
+BEAST = "2.7.2"
 RecipesBase = "1.1"
 StaticArrays = "1.2"
 Statistics = "1"
diff --git a/ext/HBEAST/HBEAST.jl b/ext/HBEAST/HBEAST.jl
new file mode 100644
index 0000000..c381b42
--- /dev/null
+++ b/ext/HBEAST/HBEAST.jl
@@ -0,0 +1,9 @@
+module HBEAST
+using BEAST
+using HMatrices
+
+include("abstractmatrix.jl")
+
+export KernelMatrix
+
+end
\ No newline at end of file
diff --git a/ext/HBEAST/abstractmatrix.jl b/ext/HBEAST/abstractmatrix.jl
new file mode 100644
index 0000000..551fa42
--- /dev/null
+++ b/ext/HBEAST/abstractmatrix.jl
@@ -0,0 +1,79 @@
+function HMatrices.KernelMatrix(
+    operator::To,
+    X::Tt,
+    Y::Ts,
+) where {To<:BEAST.IntegralOperator,Tt<:BEAST.Space,Ts<:BEAST.Space}
+    blkasm = BEAST.blockassembler(operator, X, Y)
+    function blkassembler(Z, tdata, sdata)
+        fill!(Z, 0.0)
+        @views store(v, m, n) = (Z[m, n] += v)
+        return blkasm(tdata, sdata, store)
+    end
+
+    return HMatrices.KernelMatrix{Function,typeof(X),typeof(Y),scalartype(operator)}(
+        blkassembler,
+        X,
+        Y,
+    )
+end
+
+function Base.size(K::KernelMatrix{Tf,TX,TY,T}) where {Tf,TX<:BEAST.Space,TY<:BEAST.Space,T}
+    return length(K.X), length(K.Y)
+end
+
+function Base.getindex(
+    K::KernelMatrix{Tf,TX,TY,T},
+    i::Int,
+    j::Int,
+) where {Tf,TX<:BEAST.Space,TY<:BEAST.Space,T}
+    blk = zeros(T, 1, 1)
+    K.f(blk, [i], [j])
+    return blk[1, 1]
+end
+
+function Base.getindex(
+    K::KernelMatrix{Tf,TX,TY,T},
+    i::Union{UnitRange{Int},Vector{Int}},
+    j::Union{UnitRange{Int},Vector{Int}},
+) where {Tf,TX<:BEAST.Space,TY<:BEAST.Space,T}
+    blk = zeros(T, length(i), length(j))
+    K.f(blk, i, j)
+    return blk
+end
+
+HMatrices.ClusterTree(X::BEAST.Space) = HMatrices.ClusterTree(X.pos)
+
+function HMatrices.assemble_hmatrix(
+    K::KernelMatrix{Tf,TX,TY,T};
+    rowtree = ClusterTree(K.X),
+    coltree = ClusterTree(K.Y),
+    kwargs...,
+) where {Tf,TX<:BEAST.Space,TY<:BEAST.Space,T}
+    return HMatrices.assemble_hmatrix(K, rowtree, coltree; kwargs...)
+end
+
+function HMatrices.getblock!(
+    out,
+    K::KernelMatrix{Tf,TX,TY,T},
+    irange_,
+    jrange_,
+) where {Tf,TX<:BEAST.Space,TY<:BEAST.Space,T}
+    irange = irange_ isa Colon ? axes(K, 1) : irange_
+    jrange = jrange_ isa Colon ? axes(K, 2) : jrange_
+    K.f(out, irange, jrange)
+    return out
+end
+
+function HMatrices.getblock!(
+    out,
+    K::HMatrices.PermutedMatrix{KernelMatrix{Tf,TX,TY,Tm},T},
+    irange_,
+    jrange_,
+) where {Tf,TX<:BEAST.Space,TY<:BEAST.Space,Tm,T}
+    irange = irange_ isa Colon ? axes(K, 1) : irange_
+    jrange = jrange_ isa Colon ? axes(K, 2) : jrange_
+    permuted_irange = K.rowperm[irange]
+    permuted_jrange = K.colperm[jrange]
+    K.data.f(out, permuted_irange, permuted_jrange)
+    return out
+end
diff --git a/src/HMatrices.jl b/src/HMatrices.jl
index b423be5..5e11324 100644
--- a/src/HMatrices.jl
+++ b/src/HMatrices.jl
@@ -68,6 +68,10 @@ include("triangular.jl")
 include("lu.jl")
 include("cholesky.jl")
 
+if !isdefined(Base, :get_extension) # for julia version < 1.9
+    include("../ext/HBEAST/HBEAST.jl")
+end
+
 export ClusterTree,
     CardinalitySplitter,
     DyadicSplitter,
diff --git a/test/Project.toml b/test/Project.toml
index a2203d0..c4541d3 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,6 @@
 [deps]
+BEAST = "bb4162c7-ba94-5a20-af32-d8ec4428bdd1"
+CompScienceMeshes = "3e66a162-7b8c-5da0-b8f8-124ecd2c3ae1"
 ComputationalResources = "ed09eef8-17a6-5b46-8889-db040fac31e3"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/test/extBEAST_test.jl b/test/extBEAST_test.jl
new file mode 100644
index 0000000..7456ead
--- /dev/null
+++ b/test/extBEAST_test.jl
@@ -0,0 +1,43 @@
+using BEAST
+using HMatrices
+using CompScienceMeshes
+using StaticArrays
+using LinearAlgebra
+using Test
+
+@testset "BEAST KernelMatrix" begin
+    ΓX = meshrectangle(1.0, 1.0, 0.2)
+    ΓY = translate(ΓX, SVector(3.0, 0.0, 0.0))
+
+    X = lagrangecxd0(ΓX)
+    Y = lagrangecxd0(ΓY)
+    op = Helmholtz3D.singlelayer()
+    Kref = assemble(op, X, Y)
+
+    #BEAST kernel matrix
+    K = KernelMatrix(op, X, Y)
+
+    @test Kref[1, 1] == K[1, 1]
+    @test Kref[2, 1] == K[2, 1]
+    @test Kref[1, 2] == K[1, 2]
+    @test Kref[2:5, 5:7] == K[2:5, 5:7]
+    @test size(Kref) == size(K)
+
+    splitter = HMatrices.GeometricSplitter(4)
+    rowtree = ClusterTree(X.pos, splitter)
+    coltree = ClusterTree(Y.pos, splitter)
+
+    Kperm = HMatrices.PermutedMatrix(
+        K,
+        HMatrices.loc2glob(rowtree),
+        HMatrices.loc2glob(coltree),
+    )
+
+    out = zeros(eltype(Kref), 1:size(K, 1), 1:size(K, 2))
+    HMatrices.getblock!(out, K, 1:size(K, 1), 1:size(K, 2))
+    outperm = zeros(eltype(Kref), 1:size(K, 1), 1:size(K, 2))
+    HMatrices.getblock!(outperm, Kperm, 1:size(K, 1), 1:size(K, 2))
+
+    @test out != outperm
+    @test out[Kperm.rowperm, Kperm.colperm] == outperm
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 8d3b0aa..6712f1e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -36,4 +36,7 @@ end
     include("cholesky_test.jl")
 end
 
+@safetestset "Extensions" begin
+    include("extBEAST_test.jl")
+end
 # @safetestset "DHMatrix" begin include("dhmatrix_test.jl") end

From f061ac824754f7d87ad7d96ddd370d83ad1ae8fb Mon Sep 17 00:00:00 2001
From: Joshua <joshua.tetzner@uni-rostock.de>
Date: Tue, 9 Sep 2025 18:36:35 +0200
Subject: [PATCH 2/2] Fixed formatting

---
 ext/HBEAST/HBEAST.jl  | 2 +-
 test/extBEAST_test.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/ext/HBEAST/HBEAST.jl b/ext/HBEAST/HBEAST.jl
index c381b42..fc180a9 100644
--- a/ext/HBEAST/HBEAST.jl
+++ b/ext/HBEAST/HBEAST.jl
@@ -6,4 +6,4 @@ include("abstractmatrix.jl")
 
 export KernelMatrix
 
-end
\ No newline at end of file
+end
diff --git a/test/extBEAST_test.jl b/test/extBEAST_test.jl
index 7456ead..ead485c 100644
--- a/test/extBEAST_test.jl
+++ b/test/extBEAST_test.jl
@@ -40,4 +40,4 @@ using Test
 
     @test out != outperm
     @test out[Kperm.rowperm, Kperm.colperm] == outperm
-end
\ No newline at end of file
+end