BEAST extension

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"

ext/HBEAST/HBEAST.jl

@@ -0,0 +1,9 @@
+module HBEAST
+using BEAST
+using HMatrices
+
+include("abstractmatrix.jl")
+
+export KernelMatrix
+
+end

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

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,

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"

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

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