From 25dad1f0fdef47d1f9a96bb75886baeabf5da7f0 Mon Sep 17 00:00:00 2001 From: Kyle Beggs Date: Wed, 11 Oct 2023 17:41:47 -0400 Subject: [PATCH 1/2] use safe threading --- Project.toml | 1 + src/RadialBasisFunctions.jl | 1 + src/basis/gaussian.jl | 6 +-- src/basis/inverse_multiquadric.jl | 6 +-- src/basis/polyharmonic_spline.jl | 4 +- src/interpolation.jl | 4 +- src/linalg/stencil.jl | 89 ++++++++----------------------- src/operators/operators.jl | 34 ++++-------- src/utils.jl | 6 +-- 9 files changed, 47 insertions(+), 104 deletions(-) diff --git a/Project.toml b/Project.toml index 1d39f9b..e1bbe09 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Kyle Beggs"] version = "0.1.3" [deps] +ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl index 9e9df8b..cd76d23 100644 --- a/src/RadialBasisFunctions.jl +++ b/src/RadialBasisFunctions.jl @@ -4,6 +4,7 @@ using NearestNeighbors using SymRCM using LinearAlgebra using LoopVectorization +using ChunkSplitters using SparseArrays using StaticArrays using Statistics diff --git a/src/basis/gaussian.jl b/src/basis/gaussian.jl index fdf4292..1e49565 100644 --- a/src/basis/gaussian.jl +++ b/src/basis/gaussian.jl @@ -46,11 +46,11 @@ end function Base.show(io::IO, rbf::Gaussian) print(io, "Gaussian, exp(-(ε*r)²)") - print(io, "\n└─Shape factor: ε = $(rbf.ε)") + print(io, "\n ├─Shape factor: ε = $(rbf.ε)") if rbf.poly_deg < 0 - print(io, "\n No Monomial augmentation") + print(io, "\n └─No Monomial augmentation") else - print(io, "\n└─Polynomial augmentation: degree $(rbf.poly_deg)") + print(io, "\n └─Polynomial augmentation: degree $(rbf.poly_deg)") end end diff --git a/src/basis/inverse_multiquadric.jl b/src/basis/inverse_multiquadric.jl index d3a77c5..3f774de 100644 --- a/src/basis/inverse_multiquadric.jl +++ b/src/basis/inverse_multiquadric.jl @@ -58,11 +58,11 @@ end function Base.show(io::IO, rbf::IMQ) print(io, "Inverse Multiquadrics, 1/sqrt((r*ε)²+1)") - print(io, "\n└─Shape factor: ε = $(rbf.ε)") + print(io, "\n ├─Shape factor: ε = $(rbf.ε)") if rbf.poly_deg < 0 - print(io, "\n No Monomial augmentation") + print(io, "\n └─No Monomial augmentation") else - print(io, "\n└─Polynomial augmentation: degree $(rbf.poly_deg)") + print(io, "\n └─Polynomial augmentation: degree $(rbf.poly_deg)") end end diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl index b1016ac..e293bfb 100644 --- a/src/basis/polyharmonic_spline.jl +++ b/src/basis/polyharmonic_spline.jl @@ -173,9 +173,9 @@ end function Base.show(io::IO, rbf::R) where {R<:AbstractPHS} print(io, print_basis(rbf)) if rbf.poly_deg < 0 - print(io, "\n No Monomial augmentation") + print(io, "\n └─No Monomial augmentation") else - print(io, "\n└─Polynomial augmentation: degree $(rbf.poly_deg)") + print(io, "\n └─Polynomial augmentation: degree $(rbf.poly_deg)") end end diff --git a/src/interpolation.jl b/src/interpolation.jl index a221874..027968c 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -52,8 +52,8 @@ end # pretty printing function Base.show(io::IO, op::RadialBasisInterp) println(io, "RadialBasisInterp") - println(io, " └─Input type: ", typeof(first(op.x))) - println(io, " └─Output type: ", typeof(first(op.y))) + println(io, " ├─Input type: ", typeof(first(op.x))) + println(io, " ├─Output type: ", typeof(first(op.y))) println(io, " └─Number of points: ", length(op.x)) return println( io, diff --git a/src/linalg/stencil.jl b/src/linalg/stencil.jl index 5c6cf05..ff8be3c 100644 --- a/src/linalg/stencil.jl +++ b/src/linalg/stencil.jl @@ -3,7 +3,8 @@ function _build_weightmx( data::AbstractVector{D}, centers::AbstractVector{D}, adjl::Vector{Vector{T}}, - basis::B, + basis::B; + nchunks=Threads.nthreads(), ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis} TD = eltype(first(data)) dim = length(first(data)) # dimension of data @@ -17,89 +18,45 @@ function _build_weightmx( ℒrbf = ℒ(basis) # allocate arrays to build sparse matrix - I = zeros(Int, k * length(adjl)) + Na = length(adjl) + I = zeros(Int, k * Na) J = reduce(vcat, adjl) - V = zeros(TD, k * length(adjl)) - - n_threads = Threads.nthreads() - n_threads = 1 + V = zeros(TD, k * Na) # create work arrays n = sum(sizes) - A = Symmetric[Symmetric(zeros(TD, n, n), :U) for _ in 1:n_threads] - b = Vector[zeros(TD, n) for _ in 1:n_threads] - range = Vector{UnitRange{<:Int}}(undef, n_threads) - d = Vector{Vector{eltype(data)}}(undef, n_threads) + A = Symmetric[Symmetric(zeros(TD, n, n), :U) for _ in 1:nchunks] + b = Vector[zeros(TD, n) for _ in 1:nchunks] + d = Vector{Vector{eltype(data)}}(undef, nchunks) # build stencil for each data point and store in global weight matrix - #Threads.@threads for i in eachindex(adjl) - for i in eachindex(adjl) - range[Threads.threadid()] = ((i - 1) * k + 1):(i * k) - @turbo I[range[Threads.threadid()]] .= i - d[Threads.threadid()] = data[adjl[i]] - V[range[Threads.threadid()]] = @views _build_stencil!( - A, b, Threads.threadid(), ℒrbf, ℒmon, d, centers[i], basis, mon, k - ) + Threads.@threads for (xrange, ichunk) in chunks(adjl, nchunks) + for i in xrange + I[((i - 1) * k + 1):(i * k)] .= i + d[ichunk] = data[adjl[i]] + V[((i - 1) * k + 1):(i * k)] = @views _build_stencil!( + A[ichunk], b[ichunk], ℒrbf, ℒmon, d[ichunk], centers[i], basis, mon, k + ) + end end return sparse(I, J, V, length(adjl), length(data)) end -function _build_weight_vec( - ℒ, - data::AbstractVector{D}, - centers::AbstractVector{D}, - adjl::Vector{Vector{T}}, - basis::B, -) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis} - TD = eltype(first(data)) - dim = length(first(data)) # dimension of data - nmon = binomial(dim + basis.poly_deg, basis.poly_deg) - k = length(first(adjl)) # number of data in influence/support domain - sizes = (k, nmon) - - # build monomial basis and operator - mon = MonomialBasis(dim, basis.poly_deg) - ℒmon = ℒ(mon) - ℒrbf = ℒ(basis) - - # allocate arrays to build sparse matrix - V = [zeros(TD, k) for _ in 1:length(adjl)] - - n_threads = Threads.nthreads() - - # create work arrays - n = sum(sizes) - A = Symmetric[Symmetric(zeros(TD, n, n), :U) for _ in 1:n_threads] - b = Vector[zeros(TD, n) for _ in 1:n_threads] - d = Vector{Vector{eltype(data)}}(undef, n_threads) - - # build stencil for each data point and store in global weight matrix - Threads.@threads for i in eachindex(adjl) - d[Threads.threadid()] = data[adjl[i]] - V[i] = @views _build_stencil!( - A, b, Threads.threadid(), ℒrbf, ℒmon, d, centers[i], basis, mon, k - ) - end - - return V -end - function _build_stencil!( - A::Vector{<:Symmetric}, - b::Vector{<:Vector}, - id::Int, + A::Symmetric, + b::Vector, ℒrbf, ℒmon, data::AbstractVector{D}, - center, + center::D, basis::B, mon::MonomialBasis, k::Int, ) where {D<:AbstractArray,B<:AbstractRadialBasis} - _build_collocation_matrix!(A[id], data[id], basis, mon, k) - _build_rhs!(b[id], ℒrbf, ℒmon, data[id], center, basis, k) - return (A[id] \ b[id])[1:k] + _build_collocation_matrix!(A, data, basis, mon, k) + _build_rhs!(b, ℒrbf, ℒmon, data, center, basis, k) + return (A \ b)[1:k] end function _build_collocation_matrix!( @@ -121,7 +78,7 @@ function _build_collocation_matrix!( end function _build_rhs!( - b::AbstractVector, ℒrbf, ℒmon, data::AbstractVector{D}, center, basis::B, k::K + b::AbstractVector, ℒrbf, ℒmon, data::AbstractVector{D}, center::D, basis::B, k::K ) where {D<:AbstractArray,B<:AbstractRadialBasis,K<:Int} # radial basis section @inbounds for i in eachindex(data) diff --git a/src/operators/operators.jl b/src/operators/operators.jl index 8a9bf0a..54fca7b 100644 --- a/src/operators/operators.jl +++ b/src/operators/operators.jl @@ -24,17 +24,12 @@ end # convienience constructors function RadialBasisOperator( - ℒ, - data::AbstractVector{D}, - basis::B=PHS(3; poly_deg=2); - k::T=autoselect_k(data, basis), - sparse=true, + ℒ, data::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis) ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis} - TD = eltype(D) adjl = find_neighbors(data, k) Na = length(adjl) Nd = length(data) - weights = _allocate_weights(TD, Na, Nd, k; sparse=sparse) + weights = spzeros(eltype(D), Na, Nd) return RadialBasisOperator(ℒ, weights, data, data, adjl, basis) end @@ -44,13 +39,11 @@ function RadialBasisOperator( centers::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis), - sparse=true, ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis} - TD = eltype(D) adjl = find_neighbors(data, centers, k) Na = length(adjl) Nd = length(data) - weights = _allocate_weights(TD, Na, Nd, k; sparse=sparse) + weights = spzeros(eltype(D), Na, Nd) return RadialBasisOperator(ℒ, weights, data, centers, adjl, basis) end @@ -149,19 +142,11 @@ end # update weights function _build_weights(ℒ, op::RadialBasisOperator) - if op.weights isa Vector{<:Vector} - return _build_weight_vec(ℒ, op.data, op.adjl, op.basis) - else - return _build_weightmx(ℒ, op.data, op.centers, op.adjl, op.basis) - end + return _build_weightmx(ℒ, op.data, op.centers, op.adjl, op.basis) end function _build_weights(ℒ, op::RadialBasisOperator{<:VectorValuedOperator}) - if first(op.weights) isa Vector{<:Vector} - return _build_weight_vec(ℒ, op.data, op.adjl, op.basis) - else - return _build_weightmx(ℒ, op.data, op.centers, op.adjl, op.basis) - end + return _build_weightmx(ℒ, op.data, op.centers, op.adjl, op.basis) end function update_weights!(op::RadialBasisOperator) @@ -196,11 +181,10 @@ end # pretty printing function Base.show(io::IO, op::RadialBasisOperator) println(io, "RadialBasisOperator") - println(io, " └─Operator: " * print_op(op.ℒ)) - println(io, " └─Data type: ", typeof(first(op.data))) - println(io, " └─Number of points: ", length(op.data)) - println(io, " └─Dimensions: ", length(first(op.data))) - println(io, " └─Stencil size: ", length(first(op.adjl))) + println(io, " ├─Operator: " * print_op(op.ℒ)) + println(io, " ├─Data type: ", typeof(first(op.data))) + println(io, " ├─Number of points: ", length(op.data)) + println(io, " ├─Stencil size: ", length(first(op.adjl))) return println( io, " └─Basis: ", diff --git a/src/utils.jl b/src/utils.jl index b694664..d3f2b32 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -63,7 +63,7 @@ function scale_cloud(data) return data ./ furthest_point end -_allocate_weights(m, n, k; sparse=true) = _allocate_weights(Float64, m, n, k; sparse=sparse) -function _allocate_weights(T, m, n, k; sparse=true) - return sparse ? spzeros(T, m, n) : [zeros(T, k) for _ in 1:m] +_allocate_weights(m, n, k) = _allocate_weights(Float64, m, n, k) +function _allocate_weights(T, m, n, k) + return spzeros(T, m, n) end From efc430a1cc3df8859b0d4208b123e8469c74aacb Mon Sep 17 00:00:00 2001 From: Kyle Beggs Date: Wed, 11 Oct 2023 17:46:15 -0400 Subject: [PATCH 2/2] remove sparse keyword and bump version --- Project.toml | 2 +- src/operators/gradient.jl | 16 +++++----------- 2 files changed, 6 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index e1bbe09..ed3504f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RadialBasisFunctions" uuid = "79ee0514-adf7-4479-8807-6f72ea8967e8" authors = ["Kyle Beggs"] -version = "0.1.3" +version = "0.1.4" [deps] ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e" diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl index 1a38f26..3c5121b 100644 --- a/src/operators/gradient.jl +++ b/src/operators/gradient.jl @@ -9,10 +9,7 @@ end # convienience constructors function gradient( - data::AbstractVector{D}, - basis::B=PHS(3; poly_deg=2); - k::T=autoselect_k(data, basis), - sparse=true, + data::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis) ) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int} f = ntuple(length(first(data))) do dim return let dim = dim @@ -20,7 +17,7 @@ function gradient( end end ℒ = Gradient(f) - return RadialBasisOperator(ℒ, data, basis; k=k, sparse=sparse) + return RadialBasisOperator(ℒ, data, basis; k=k) end function gradient( @@ -28,7 +25,6 @@ function gradient( centers::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis), - sparse=true, ) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int} f = ntuple(length(first(data))) do dim return let dim = dim @@ -36,7 +32,7 @@ function gradient( end end ℒ = Gradient(f) - return RadialBasisOperator(ℒ, data, centers, basis; k=k, sparse=sparse) + return RadialBasisOperator(ℒ, data, centers, basis; k=k) end function RadialBasisOperator( @@ -44,12 +40,11 @@ function RadialBasisOperator( data::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis), - sparse=true, ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis} TD = eltype(D) adjl = find_neighbors(data, k) N = length(adjl) - weights = ntuple(_ -> _allocate_weights(TD, N, N, k; sparse=sparse), length(ℒ.ℒ)) + weights = ntuple(_ -> _allocate_weights(TD, N, N, k), length(ℒ.ℒ)) return RadialBasisOperator(ℒ, weights, data, data, adjl, basis) end @@ -59,13 +54,12 @@ function RadialBasisOperator( centers::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis), - sparse=true, ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis} TD = eltype(D) adjl = find_neighbors(data, centers, k) Na = length(adjl) Nd = length(data) - weights = ntuple(_ -> _allocate_weights(TD, Na, Nd, k; sparse=sparse), length(ℒ.ℒ)) + weights = ntuple(_ -> _allocate_weights(TD, Na, Nd, k), length(ℒ.ℒ)) return RadialBasisOperator(ℒ, weights, data, centers, adjl, basis) end