Merge pull request #15 from kylebeggs/feature/multithreading

Use safe multhreading

Project.toml

@@ -1,9 +1,10 @@
 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"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/RadialBasisFunctions.jl

@@ -4,6 +4,7 @@ using NearestNeighbors
 using SymRCM
 using LinearAlgebra
 using LoopVectorization
+using ChunkSplitters
 using SparseArrays
 using StaticArrays
 using Statistics

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
 

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
 

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
 

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,

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)

src/operators/gradient.jl

@@ -9,47 +9,42 @@ 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
  x -> ∂(x, 1, dim)
  end
  end
  ℒ = Gradient(f)
- return RadialBasisOperator(ℒ, data, basis; k=k, sparse=sparse)
+ return RadialBasisOperator(ℒ, data, basis; k=k)
 end
 
 function gradient(
  data::AbstractVector{D},
  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
  x -> ∂(x, 1, dim)
  end
  end
  ℒ = Gradient(f)
- return RadialBasisOperator(ℒ, data, centers, basis; k=k, sparse=sparse)
+ return RadialBasisOperator(ℒ, data, centers, basis; k=k)
 end
 
 function RadialBasisOperator(
  ℒ::Gradient,
  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
 

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: ",

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