API Reference

Full docstrings for every exported name, plus internals for developers extending the package.

Macros

RadialBasisFunctions.@operator Macro

@operator expr

Create an operator from mathematical notation. Returns an AbstractOperator that can be called directly with data points to build a RadialBasisOperator.

Recognized symbols

• ∇² / Δ — Laplacian

• ∂(dim) — first partial derivative in dimension dim

• ∂²(dim) — second partial derivative in dimension dim

• ∂(dim1, dim2) — mixed partial derivative ∂²f/(∂xᵢ ∂xⱼ)

• ∇ ⋅ (κ * ∇) — anisotropic diffusion (scalar or vector κ)

• c ⋅ ∇ — advection operator (vector c)

• f / I — Identity operator

• Everything else — scalar coefficient

Examples

helm = (@operator ∇² + k^2 * f)(x)
aniso = (@operator κx * ∂²(1) + κy * ∂²(2))(x)
advdiff = (@operator ν * ∇² - c ⋅ ∇)(x)
diff = (@operator ∇ ⋅ (κ * ∇))(x)

source

Exported Functions

RadialBasisFunctions.PHS Method

function PHS(n::T=3; poly_deg::T=2) where {T<:Int}

Convenience constructor for polyharmonic splines.

Arguments

• n: Order of the spline (1, 3, 5, or 7). Higher = smoother.

• poly_deg: Polynomial augmentation degree (default: 2 for quadratic).

See also: IMQ, Gaussian

source

RadialBasisFunctions.autoselect_k Method

autoselect_k(data::Vector, basis<:AbstractRadialBasis)

Return the stencil size k (number of nearest neighbors per stencil) used by default when constructing operators: min(N, max(2*binomial(m+d, d), 2d+1)) where N = length(data), m is the polynomial degree of basis, and d is the spatial dimension.

See Bayona, 2017 - https://doi.org/10.1016/j.jcp.2016.12.008

source

RadialBasisFunctions.classify_stencil Method

classify_stencil(is_boundary, boundary_conditions, eval_idx,
                neighbors, global_to_boundary)

Classify stencil type for dispatch in kernel execution.

source

RadialBasisFunctions.curl Method

curl(data, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a curl operator and applies it to vector field x.

For repeated evaluations on the same points, prefer creating the operator once with curl(data) and calling it via functor syntax op(x).

source

RadialBasisFunctions.curl Method

curl(data; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for the curl (∇×u).

Only defined for 2D and 3D data.

• 2D: Returns scalar field ∂u₂/∂x₁ − ∂u₁/∂x₂

• 3D: Returns vector field with standard curl components

Arguments

• data: Vector of data points

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

# 2D curl (scalar output)
points = rand(SVector{2,Float64}, 1000)
curl_op = curl(points)
u = hcat(-getindex.(points, 2), getindex.(points, 1))  # u = (-y, x)
ω = curl_op(u)  # ≈ 2.0 everywhere

# 3D curl (vector output)
points3d = rand(SVector{3,Float64}, 1000)
curl_op3d = curl(points3d)
u3d = hcat(-getindex.(points3d, 2), getindex.(points3d, 1), zeros(1000))
ω3d = curl_op3d(u3d)  # ≈ [0, 0, 2] everywhere

See also: divergence, gradient, jacobian

source

RadialBasisFunctions.custom Method

custom(data, ℒ::Function; rank=<auto>, basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)
custom(data, op::AbstractOperator; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator with a custom operator.

Arguments

• data: Vector of data points

• ℒ: Custom function that accepts a basis and returns a callable (x, xᵢ) -> value

• op: An AbstractOperator (e.g. from @operator or operator algebra)

Keyword Arguments

• rank::Int: Tensor rank added to the output (0 = rank-preserving, 1 = rank+1). Auto-inferred when omitted: from the type parameter for AbstractOperator, or by probing the closure for Function (tuple return → rank 1, scalar → rank 0).

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

# Using @operator macro — call the operator directly with data points
op = (@operator ∇² + k² * f)(data)

# Custom operator that returns the basis function itself (rank-preserving)
op = custom(data, basis -> (x, xᵢ) -> basis(x, xᵢ))

# Custom second partial derivative ∂²f/∂x₁² using the ∂² functor
op = custom(data, basis -> ∂²(basis, 1))

source

RadialBasisFunctions.derivative_order Method

derivative_order(op) -> Union{Int, Missing}

Total order of differentiation. Useful for estimating required polynomial degree and stencil size. Returns missing for Custom operators, whose derivative order cannot be inferred from the wrapped function.

source

RadialBasisFunctions.directional Method

directional(data, v; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for the directional derivative (∇f⋅v).

Arguments

• data: Vector of data points

• v: Direction vector. Can be:

  • A single vector of length Dim (constant direction)

  • A vector of vectors matching length(data) (spatially-varying direction)

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

# Constant direction
∂_x = directional(data, [1.0, 0.0])

# Spatially-varying direction (e.g., radial)
normals = normalize.(data)
∂_n = directional(data, normals)

See also: gradient, partial, laplacian

source

RadialBasisFunctions.divergence Method

divergence(data, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a divergence operator and applies it to vector field x.

For repeated evaluations on the same points, prefer creating the operator once with divergence(data) and calling it via functor syntax op(x).

source

RadialBasisFunctions.divergence Method

divergence(data; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for the divergence (∇⋅u).

Arguments

• data: Vector of data points

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

points = rand(SVector{2,Float64}, 1000)
div_op = divergence(points)

# Vector field as matrix (N × D)
u = hcat(sin.(getindex.(points, 1)), cos.(getindex.(points, 2)))
div_u = div_op(u)  # Vector (N,)

See also: curl, gradient, jacobian

source

RadialBasisFunctions.gradient Method

gradient(data, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a gradient operator and applies it to scalar field x.

For repeated evaluations, prefer creating the operator once with gradient(data).

Examples

points = rand(SVector{2,Float64}, 1000)
u = sin.(getindex.(points, 1))
∇u = gradient(points, u)  # One-shot gradient computation

source

RadialBasisFunctions.gradient Method

gradient(data; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for computing gradients of scalar fields.

This is a convenience alias for jacobian. The gradient of a scalar field is mathematically the Jacobian (a 1×D row vector, returned as a length-D vector per point).

Arguments

• data: Vector of points (e.g., Vector{SVector{2,Float64}})

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

points = rand(SVector{2,Float64}, 1000)
op = gradient(points)

u = sin.(getindex.(points, 1))
∇u = op(u)  # Matrix (1000 × 2)
∂u_∂x = ∇u[:, 1]
∂u_∂y = ∇u[:, 2]

See also: jacobian

source

RadialBasisFunctions.hessian Method

hessian(data, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a Hessian operator and applies it to field x.

source

RadialBasisFunctions.hessian Method

hessian(data; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for computing the full Hessian matrix.

Arguments

• data: Vector of points (e.g., Vector{SVector{2,Float64}})

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

points = rand(SVector{2,Float64}, 1000)
op = hessian(points)

# Scalar field → Hessian tensor
u = sin.(getindex.(points, 1)) .* cos.(getindex.(points, 2))
H = op(u)  # Array (1000 × 2 × 2)
# H[:, 1, 1] = ∂²u/∂x², H[:, 1, 2] = ∂²u/∂x∂y, etc.

See also: jacobian, laplacian, mixed_partial

source

RadialBasisFunctions.is_antisymmetric Method

is_antisymmetric(op) -> Bool

Whether the operator produces an anti-symmetric output tensor (Aᵢⱼ = −Aⱼᵢ).

source

RadialBasisFunctions.is_self_adjoint Method

is_self_adjoint(op) -> Bool

Whether the operator is self-adjoint (⟨ℒu, v⟩ = ⟨u, ℒv⟩). Self-adjoint operators produce symmetric weight matrices, which matters for solver selection and eigenvalue problems.

source

RadialBasisFunctions.is_symmetric Method

is_symmetric(op) -> Bool

Whether the operator produces a symmetric output tensor. Symmetric operators can exploit storage savings (e.g., only storing upper-triangular entries).

source

RadialBasisFunctions.jacobian Method

jacobian(data, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a Jacobian operator and applies it to field x.

For repeated evaluations on the same points, prefer creating the operator once with jacobian(data) and calling it via functor syntax op(x).

Arguments

• data: Vector of points

• x: Field values to differentiate

Keyword Arguments

• basis: RBF basis function (default: PHS(3; poly_deg=2))

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

Examples

points = rand(SVector{2,Float64}, 1000)
u = sin.(getindex.(points, 1))
∇u = jacobian(points, u)  # One-shot gradient computation

source

RadialBasisFunctions.jacobian Method

jacobian(data; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for computing Jacobians (or gradients for scalar fields).

The Jacobian is the fundamental differential operator. For a scalar field, it computes the gradient. For a vector field, it computes the full Jacobian matrix. The spatial dimension is automatically inferred from the data.

Arguments

• data: Vector of points (e.g., Vector{SVector{2,Float64}})

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

points = rand(SVector{2,Float64}, 1000)
op = jacobian(points)

# Scalar field → gradient
u = sin.(getindex.(points, 1))
∇u = op(u)  # Matrix (1000 × 2)

# Vector field → Jacobian matrix
v = hcat(u, cos.(getindex.(points, 2)))
J = op(v)  # Array (1000 × 2 × 2)

See also: gradient

source

RadialBasisFunctions.laplacian Method

laplacian(data; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for the Laplacian operator (∇²f).

Arguments

• data: Vector of data points

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

# Basic usage
op = laplacian(data)

# With custom basis
op = laplacian(data; basis=PHS(5; poly_deg=3))

# With different evaluation points
op = laplacian(data; eval_points=eval_pts)

See also: partial, gradient, directional

source

RadialBasisFunctions.mixed_partial Method

mixed_partial(data, dim1, dim2, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a mixed partial operator and applies it to field x.

source

RadialBasisFunctions.mixed_partial Method

mixed_partial(data, dim1, dim2; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for the mixed partial derivative ∂²f/(∂xᵢ ∂xⱼ).

Arguments

• data: Vector of data points

• dim1: First dimension index

• dim2: Second dimension index

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

∂²xy = mixed_partial(data, 1, 2)

See also: partial, hessian

source

RadialBasisFunctions.normal_derivative Method

normal_derivative(data, normals, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a normal derivative operator and applies it to field x.

For repeated evaluations on the same points, prefer creating the operator once with normal_derivative(data, normals) and calling it via functor syntax op(x).

source

RadialBasisFunctions.normal_derivative Method

normal_derivative(data, normals; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for the normal derivative (∇f⋅n̂).

The input normals are automatically normalized to unit vectors. This is a convenience wrapper around directional.

Arguments

• data: Vector of data points

• normals: Normal vectors at each point (will be normalized)

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

points = rand(SVector{2,Float64}, 1000)
normals = normalize.(points)  # radial normals
∂ₙ = normal_derivative(points, normals)
result = ∂ₙ(sin.(getindex.(points, 1)))

See also: directional, gradient

source

RadialBasisFunctions.output_rank Method

output_rank(op::AbstractOperator{N}) -> Int

Tensor rank added to the output by this operator. Returns N from the type parameter.

source

RadialBasisFunctions.partial Method

partial(data, order, dim; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for a partial derivative.

Arguments

• data: Vector of data points

• order: Derivative order (1 or 2)

• dim: Dimension index to differentiate

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

# First derivative in x-direction
∂x = partial(data, 1, 1)

# Second derivative in y-direction
∂²y = partial(data, 2, 2; basis=PHS(5; poly_deg=4))

See also: laplacian, gradient, directional

source

RadialBasisFunctions.regrid Method

regrid(data, eval_points; basis=PHS(3; poly_deg=2), k, adjl)

Build a RadialBasisOperator for interpolating from data points to eval_points.

Arguments

• data: Source data points

• eval_points: Target evaluation points

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

Examples

# Interpolate from coarse grid to fine grid
coarse = rand(SVector{2,Float64}, 100)
fine = rand(SVector{2,Float64}, 1000)
op = regrid(coarse, fine)

# Apply to field values
u_coarse = sin.(getindex.(coarse, 1))
u_fine = op(u_coarse)

See also: Interpolator

source

RadialBasisFunctions.rotation_rate Method

rotation_rate(data, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a rotation rate operator and applies it to vector field x.

source

RadialBasisFunctions.rotation_rate Method

rotation_rate(data; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for the anti-symmetric rotation rate tensor ωᵢⱼ = ½(∂uᵢ/∂xⱼ − ∂uⱼ/∂xᵢ).

Arguments

• data: Vector of data points

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

points = rand(SVector{2,Float64}, 1000)
ω_op = rotation_rate(points)

# Solid body rotation: u = (-y, x) → ω₁₂ = -1
u = hcat(-getindex.(points, 2), getindex.(points, 1))
ω = ω_op(u)  # Array (1000 × 2 × 2), ω₁₂ ≈ -1.0

See also: strain_rate, jacobian, curl

source

RadialBasisFunctions.strain_rate Method

strain_rate(data, x; basis=PHS(3; poly_deg=2), k, adjl)

One-shot convenience function that creates a strain rate operator and applies it to vector field x.

source

RadialBasisFunctions.strain_rate Method

strain_rate(data; basis=PHS(3; poly_deg=2), eval_points=data, k, adjl, hermite)

Build a RadialBasisOperator for the symmetric strain rate tensor εᵢⱼ = ½(∂uᵢ/∂xⱼ + ∂uⱼ/∂xᵢ).

Arguments

• data: Vector of data points

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

Examples

points = rand(SVector{2,Float64}, 1000)
ε_op = strain_rate(points)

# Vector field as matrix (N × D)
u = hcat(getindex.(points, 2), getindex.(points, 1))  # u = (y, x)
ε = ε_op(u)  # Array (1000 × 2 × 2), ε₁₂ = ε₂₁ ≈ 1.0

See also: rotation_rate, jacobian, divergence

source

RadialBasisFunctions.update_hermite_stencil_data! Method

update_hermite_stencil_data!(hermite_data, global_data, neighbors,
                             is_boundary, boundary_conditions, normals,
                             global_to_boundary[, eval_point])

Populate local Hermite stencil data from global arrays. Used within kernels to extract boundary info for specific neighborhoods.

When eval_point is given, caches its local index (first stencil point equal by value, 0 if absent) in hermite_data.eval_local_idx so RHS assembly avoids a per-operator search.

source

RadialBasisFunctions.update_weights! Method

update_weights!(op::RadialBasisOperator)

Rebuild the stencil weights of op from its current data and eval_points, then mark the weight cache as valid. For StencilWeights storage this is a pure in-place write of the value matrix — no sparse reassembly.

source

RadialBasisFunctions.weights Method

weights(op::RadialBasisOperator)

Return the stencil weights of op — a StencilWeights (dense stencil-wise ELL storage) for scalar-valued (rank-0) operators, or a tuple of StencilWeights for gradient-family operators. Rebuilds the weights first if the cache is stale.

This is the supported accessor over reaching into the weights field directly. For a SparseMatrixCSC (global system assembly, external sparse libraries), use sparse(op).

source

RadialBasisFunctions.∂virtual Method

function ∂virtual(data, eval_points, dim, Δ, basis; backward=false, k=autoselect_k(data, basis))

Builds a virtual RadialBasisOperator which will be evaluated at eval_points where the operator is the partial derivative with respect to dim. Virtual operators interpolate the data to structured points at a distance Δ for which standard finite difference formulas can be applied. The result is a standard RadialBasisOperator, so weight caching and update_weights! come for free.

Defaults to a forward difference (backward=false). Note the convenience form without eval_points instead defaults to a backward difference (backward=true).

source

RadialBasisFunctions.∂virtual Method

function ∂virtual(data, dim, Δ, basis; backward=true, k=autoselect_k(data, basis))

Builds a virtual RadialBasisOperator which will be evaluated at the input points (data) where the operator is the partial derivative with respect to dim. Virtual operators interpolate the data to structured points at a distance Δ for which standard finite difference formulas can be applied. The result is a standard RadialBasisOperator, so weight caching and update_weights! come for free.

Defaults to a backward difference (backward=true), unlike the form with explicit eval_points, which defaults to forward (backward=false).

source

RadialBasisFunctions.AbstractBasis Type

abstract type AbstractBasis end

source

RadialBasisFunctions.AbstractOperator Type

AbstractOperator{N}

Abstract supertype for differential operators.

The parameter N is the tensor rank added to the output:

• N=0: rank-preserving (e.g., Partial, Laplacian, Directional)

• N=1: adds one trailing dimension (e.g., Jacobian)

• N=2: adds two trailing dimensions (e.g., Hessian)

Subtypes must implement (op::MyOp)(basis::AbstractBasis) to return a callable applied to the basis.

All AbstractOperator subtypes are callable with data points to construct a RadialBasisOperator:

op = Laplacian()(points)
op = (@operator ∇² + k² * f)(points; basis=PHS(5; poly_deg=3))

source

RadialBasisFunctions.AbstractPHS Type

abstract type AbstractPHS <: AbstractRadialBasis

Supertype of all Polyharmonic Splines.

source

RadialBasisFunctions.AbstractRadialBasis Type

abstract type AbstractRadialBasis <: AbstractBasis end

source

RadialBasisFunctions.BoundaryCondition Type

BoundaryCondition{T}

Unified boundary condition representation: Bu = α_u + β_∂ₙu

Special cases:

• Dirichlet: α=1, β=0

• Neumann: α=0, β=1

• Robin: α≠0, β≠0

• Internal: α=0, β=0 (sentinel for interior points)

source

RadialBasisFunctions.Curl Type

Curl{Dim} <: AbstractJacobianOperator{Dim,0}

Operator for the curl of a vector field (∇×u).

• 2D: Matrix (N×2) → Vector (N): ∂u₂/∂x₁ − ∂u₁/∂x₂

• 3D: Matrix (N×3) → Matrix (N×3): standard curl vector

Only defined for Dim ∈ {2, 3}. Weights are stored as NTuple{Dim, StencilWeights}, reusing the Jacobian weight-building infrastructure.

source

RadialBasisFunctions.Custom Type

Custom{N, F<:Function} <: AbstractOperator{N}

Custom operator that applies a user-defined function to basis functions. The function ℒ should accept a basis and return a callable (x, xᵢ) -> value.

N is the tensor rank added to the output: Custom{0}(ℒ) for rank-preserving or Custom{1}(ℒ) for rank+1.

source

RadialBasisFunctions.Directional Type

Directional{Dim,T} <: AbstractOperator{0}

Operator for the directional derivative (∇f⋅v), the inner product of the gradient and a direction vector.

source

RadialBasisFunctions.Divergence Type

Divergence{Dim} <: AbstractJacobianOperator{Dim,0}

Operator for the divergence of a vector field (∇⋅u = ∑ᵢ ∂uᵢ/∂xᵢ).

Takes a vector field (Matrix N×D) as input and produces a scalar field (Vector N). Weights are stored as NTuple{Dim, StencilWeights}, one partial derivative matrix per spatial dimension, reusing the Jacobian weight-building infrastructure.

source

RadialBasisFunctions.Gaussian Type

Gaussian(ε=1; poly_deg=2)

Gaussian radial basis function: ϕ(r) = exp(-(εr)²)

Arguments

• ε: Shape parameter (must be > 0). Smaller values = wider basis.

• poly_deg: Polynomial augmentation degree (default: 2).

See also: PHS, IMQ

source

RadialBasisFunctions.HermiteStencilData Type

HermiteStencilData{T}

Local stencil data for Hermite interpolation with boundary conditions.

Fields:

• data: Coordinates of k stencil points

• is_boundary: Boolean flags for each point

• boundary_conditions: BC for each point (use Internal() for interior)

• normals: Normal vectors (zero for interior points)

• poly_workspace: Pre-allocated buffer for polynomial operations (avoids allocations in hot path)

• normal_workspace: Pre-allocated scratch for in-place normal-derivative evaluation

• eval_local_idx: Local index of the evaluation point among the stencil points (set by update_hermite_stencil_data!; 0 when absent/unset)

Note: For interior points (is_boundary[i] == false), boundary_conditions[i] and normals[i] contain sentinel values and should not be accessed.

source

RadialBasisFunctions.HermiteStencilData Method

Pre-allocation constructor for HermiteStencilData

source

RadialBasisFunctions.Hessian Type

Hessian{Dim} <: AbstractOperator{2}

Operator for computing the full Hessian matrix ∂²f/(∂xᵢ ∂xⱼ) at each point.

The Hessian is a rank-2 operator: it adds two trailing dimensions of size D (the spatial dimension) to the output. Only D*(D+1)/2 unique weight matrices are stored (upper-triangular entries); symmetric entries are filled via copyto!.

Input/Output Shapes

• Scalar field Vector{T} (N,) → Hessian Array{T,3} (N_eval × D × D)

• Vector field Matrix{T} (N × D_in) → Array{T,4} (N_eval × D_in × D × D)

source

RadialBasisFunctions.IMQ Type

IMQ(ε=1; poly_deg=2)

Inverse Multiquadric radial basis function: ϕ(r) = 1/√((εr)² + 1)

Arguments

• ε: Shape parameter (must be > 0). Smaller values = flatter basis.

• poly_deg: Polynomial augmentation degree (default: 2).

See also: PHS, Gaussian

source

RadialBasisFunctions.Identity Type

Identity <: AbstractOperator{0}

Identity operator — returns the basis function unchanged. Useful in operator algebra to represent the function itself (e.g., Laplacian() + k² * Identity() for Helmholtz).

source

RadialBasisFunctions.Interpolator Type

struct Interpolator

Construct a radial basis interpolation.

source

RadialBasisFunctions.Interpolator Method

function Interpolator(x, y, basis::B=PHS())

Construct a radial basis interpolator.

See also: regrid for local stencil-based interpolation between point sets.

source

RadialBasisFunctions.Jacobian Type

Jacobian{Dim} <: AbstractJacobianOperator{Dim,1}

Operator type for computing Jacobians (and gradients as a special case).

The Jacobian is the fundamental differential operator that computes all partial derivatives. When applied to a scalar field, it produces the gradient. When applied to a vector field, it produces the full Jacobian matrix.

Differentiation increases tensor rank by 1. The output gains a trailing dimension of size D (the spatial dimension).

Input/Output Shapes

• Scalar field Vector{T} (N,) → Gradient Matrix{T} (N_eval × D)

• Vector field Matrix{T} (N × D) → Jacobian Array{T,3} (N_eval × D × D)

• Matrix field Array{T,3} (N × D × D) → 3-tensor Array{T,4} (N_eval × D × D × D)

• General: input shape (N, dims...) → output shape (N_eval, dims..., D)

source

RadialBasisFunctions.Laplacian Type

Laplacian <: AbstractOperator{0}

Operator for the sum of the second derivatives w.r.t. each independent variable (∇²f).

source

RadialBasisFunctions.MixedPartial Type

MixedPartial{T<:Int} <: AbstractOperator{0}

Operator for the mixed second partial derivative ∂²f/(∂xᵢ ∂xⱼ). When dim1 == dim2, delegates to Partial(2, dim).

source

RadialBasisFunctions.MonomialBasis Type

struct MonomialBasis{Dim,Deg} <: AbstractBasis

Dim dimensional monomial basis of order Deg.

source

RadialBasisFunctions.PHS1 Type

struct PHS1{T<:Int} <: AbstractPHS

Polyharmonic spline radial basis function: 

source

RadialBasisFunctions.PHS3 Type

struct PHS3{T<:Int} <: AbstractPHS

Polyharmonic spline radial basis function: 

source

RadialBasisFunctions.PHS5 Type

struct PHS5{T<:Int} <: AbstractPHS

Polyharmonic spline radial basis function: 

source

RadialBasisFunctions.PHS7 Type

struct PHS7{T<:Int} <: AbstractPHS

Polyharmonic spline radial basis function: 

source

RadialBasisFunctions.Partial Type

Partial{T<:Int} <: AbstractOperator{0}

Operator for a partial derivative of specified order with respect to a dimension.

source

RadialBasisFunctions.RadialBasisOperator Type

struct RadialBasisOperator

Operator of data using a radial basis with potential monomial augmentation.

source

RadialBasisFunctions.RadialBasisOperator Method

RadialBasisOperator(ℒ, data; eval_points, basis, k, adjl, hermite, device)

Unified constructor with keyword arguments.

Arguments

• ℒ: The operator type (e.g., Laplacian(), Partial(1, 2))

• data: Vector of data points

Keyword Arguments

• basis: RBF basis (default: PHS(3; poly_deg=2))

• eval_points: Evaluation points (default: data)

• k: Stencil size (default: autoselect_k(data, basis))

• adjl: Adjacency list (default: computed via find_neighbors)

• hermite: Optional NamedTuple for Hermite interpolation

• device: KernelAbstractions backend for weight computation (default: auto-detected from data via get_backend)

The hermite NamedTuple has fields is_boundary::AbstractVector{Bool} (a BitVector works), bc::AbstractVector{<:BoundaryCondition}, and normals::AbstractVector{<:AbstractVector}.

Examples

# Preferred: call operator directly with data points
op = Laplacian()(data)
op = Laplacian()(data; basis=PHS(5; poly_deg=3), k=40)
op = (@operator ∇² + k² * f)(data)

# Explicit constructor form
op = RadialBasisOperator(Laplacian(), data)

# With different evaluation points
op = Laplacian()(data; eval_points=eval_pts)

# With Hermite boundary conditions
op = RadialBasisOperator(Laplacian(), data;
    hermite=(is_boundary=is_bound, bc=boundary_conds, normals=normal_vecs))

# With explicit device
using KernelAbstractions
op = RadialBasisOperator(Laplacian(), data; device=CPU())

source

RadialBasisFunctions.RadialBasisOperator Method

function (op::RadialBasisOperator)(y, x)

Evaluate the operator at x in-place and store the result in y.

source

RadialBasisFunctions.RadialBasisOperator Method

function (op::RadialBasisOperator)(x)

Evaluate the operator at x.

source

RadialBasisFunctions.Regrid Type

Regrid

Operator for interpolating from one set of points to another.

source

RadialBasisFunctions.RotationRate Type

RotationRate{Dim} <: AbstractJacobianOperator{Dim,0}

Operator for the anti-symmetric rotation rate tensor ωᵢⱼ = ½(∂uᵢ/∂xⱼ − ∂uⱼ/∂xᵢ).

Takes a vector field (Matrix N×D) as input and produces an anti-symmetric tensor (Array N_eval×D×D). Diagonal entries are zero. Weights are stored as NTuple{Dim, StencilWeights}, reusing the Jacobian weight-building infrastructure.

source

RadialBasisFunctions.ScaledOperator Type

ScaledOperator{N, T<:Number, O<:AbstractOperator{N}} <: AbstractOperator{N}

An operator multiplied by a scalar coefficient. Created via α * op or op * α.

source

RadialBasisFunctions.StencilWeights Type

StencilWeights{T, E <: EllSparse.SellMatrix{T}} <: AbstractMatrix{T}

ELL-format stencil weight storage for RadialBasisOperators — a thin policy wrapper over EllSparse.SellMatrix, which owns the storage layout and the CPU/GPU apply kernels. The wrapper keeps the RBF-specific policy: sparse-mixing algebra that never densifies, the \ guardrail through sparse, the frozen-structure copyto! refusal, and the parent contract.

parent(W)[l, i] is the l-th stencil weight of eval point i (column order matches the adjacency list adjl) and the frozen index structure records the data point it multiplies. The logical size is (N_eval, N_data): W * x computes y[i] = Σₗ vals[l, i] * x[idx[l, i]] with one dense length-k dot product per eval point — row-parallel on CPU threads or GPU.

The index structure is frozen after construction (update_weights! rewrites only the values), and all components of a gradient-family operator share one EllSparse.SellStructure — one index matrix and one precomputed transpose map — which makes the adjoint apply W' * x a deterministic row-parallel gather. Dirichlet boundary rows are padded to k slots: weight 1 at slot 1 with index eval_idx, zeros (with the same in-bounds index) elsewhere.

Use sparse / SparseMatrixCSC to convert to a sparse matrix for global system assembly; duplicate padded indices combine, so Dirichlet rows convert to 1-entry identity rows exactly as the previous sparse storage produced.

source

RadialBasisFunctions.StrainRate Type

StrainRate{Dim} <: AbstractJacobianOperator{Dim,0}

Operator for the symmetric strain rate tensor εᵢⱼ = ½(∂uᵢ/∂xⱼ + ∂uⱼ/∂xᵢ).

Takes a vector field (Matrix N×D) as input and produces a symmetric tensor (Array N_eval×D×D). Weights are stored as NTuple{Dim, StencilWeights}, reusing the Jacobian weight-building infrastructure.

source

RadialBasisFunctions.SumOperator Type

SumOperator{N, Ops<:Tuple{Vararg{AbstractOperator}}} <: AbstractOperator{N}

Sum of operators sharing output rank N. Created via op1 + op2 or op1 - op2 (subtraction stores -1 * op2); nested sums flatten into a single n-ary term tuple.

source

RadialBasisFunctions.VirtualPartial Type

VirtualPartial{T<:Real} <: AbstractOperator{0}

Operator for the first partial derivative with respect to dim, computed virtually: the data is interpolated (regridded) to points shifted by Δ along dim and a one-sided finite difference is taken between the shifted and unshifted interpolants.

source

Private

Base.:* Method

α * op::RadialBasisOperator
op::RadialBasisOperator * α
op::RadialBasisOperator / α
-op::RadialBasisOperator

Scale a built operator by a scalar. The weights are scaled directly (no re-collocation) and the symbolic operator is wrapped in a ScaledOperator, so α * laplacian(x) is equivalent to (@operator α * ∇²)(x) up to floating-point associativity.

source

Base.eltype Method

Base.eltype(op::RadialBasisOperator)

Return the element type of the operator's weight matrix.

source

Base.parent Method

parent(W::StencilWeights)

For host-resident weights — the only kind construction and AD produce — return the dense k × N_eval stencil-major weight value matrix backing W (zero-copy): the supported handle for reading or mutating weight values in place and for AD losses over built weights (e.g. sum(parent(W) .^ 2)). Device-adapted weights return the device-layout values array — treat it as opaque storage and go through copyto!, sparse, or Matrix instead.

source

Base.size Method

Base.size(op::RadialBasisOperator{<:AbstractOperator{0}})

Return the dimensions (m, n) of the operator's weight matrix. Only defined for rank-0 operators with scalar weights.

source

LinearAlgebra.mul! Method

LinearAlgebra.mul!(y, op::RadialBasisOperator, x)
LinearAlgebra.mul!(y, op::RadialBasisOperator, x, α, β)

Standard linear algebra interface for operator application: y = op * x or y = α * op * x + β * y.

This makes RadialBasisOperator composable with iterative solvers (Krylov.jl), LinearMaps.jl, and any generic Julia code expecting a mul!-compatible linear map.

Only defined for rank-0 operators with scalar (non-tuple) weights, since mul! expects a single linear map y = A * x.

source

RadialBasisFunctions._any_boundary Method

Early-exit scan; avoids materializing is_boundary[neighbors] per stencil.

source

RadialBasisFunctions._assemble_weights Method

_assemble_weights(vals_list, idx, N_data, num_ops)

Wrap the filled ELL buffers into StencilWeights — one per operator component, all sharing the same neighbor-index matrix.

source

RadialBasisFunctions._backward_partial_poly_1d! Method

Backward pass for polynomial section in 1D.

source

RadialBasisFunctions._backward_partial_poly_2d! Method

Backward pass for polynomial section in 2D.

source

RadialBasisFunctions._backward_partial_poly_3d! Method

Backward pass for polynomial section in 3D.

source

RadialBasisFunctions._backward_partial_polynomial_section! Method

_backward_partial_polynomial_section!(Δeval_point, Δb, k, nmon, dim, num_ops)

Backward pass through the polynomial section of the RHS for Partial operator.

For monomials in 2D with poly_deg=2 (1, x, y, xy, x², y²): ∂/∂x gives: 0, 1, 0, y, 2x, 0

The gradients of these w.r.t. eval_point are: ∂(0)/∂(x,y) = (0, 0) ∂(1)/∂(x,y) = (0, 0) ∂(0)/∂(x,y) = (0, 0) ∂(y)/∂(x,y) = (0, 1) -> b[4] contributes to Δeval_point[2] ∂(2x)/∂(x,y) = (2, 0) -> b[5] contributes 2 to Δeval_point[1] ∂(0)/∂(x,y) = (0, 0)

This is equivalent to computing the mixed second derivatives ∂²pⱼ/∂x[dim]∂x[d].

source

RadialBasisFunctions._build_collocation_matrix! Method

_build_collocation_matrix!(A, data, basis, mon, k)

Build RBF collocation matrix. Works for both interior stencils (AbstractVector data) and Hermite stencils (HermiteStencilData) via dispatch helpers.

Matrix structure:

┌─────────────────┬─────────┐
│  Φ(xᵢ, xⱼ)      │ P(xᵢ)   │  k×k RBF + k×nmon polynomial
├─────────────────┼─────────┤
│  P(xⱼ)ᵀ         │   0     │  nmon×k poly + nmon×nmon zero
└─────────────────┴─────────┘

For Hermite stencils with Neumann/Robin conditions, basis functions are modified via _rbf_entry and _poly_entry! dispatch to maintain matrix symmetry.

source

RadialBasisFunctions._build_rhs! Method

_build_rhs!(b, ℒrbf::Tuple, ℒmon::Tuple, data, eval_point, basis, mon, k)

Build RHS vector for multiple operators. Works for both interior stencils (AbstractVector) and Hermite stencils (HermiteStencilData) via dispatch helpers.

source

RadialBasisFunctions._build_rhs! Method

_build_rhs!(b, ℒrbf, ℒmon, data, eval_point, basis, mon, k)

Build RHS vector for single operator. Works for both interior stencils (AbstractVector) and Hermite stencils (HermiteStencilData) via dispatch helpers.

source

RadialBasisFunctions._build_stencil! Method

_build_stencil!(λ, A, b, ℒrbf, ℒmon, data, eval_point, basis, mon, k)

Assemble complete stencil: build collocation matrix, build RHS, solve for weights in the pre-allocated λ buffer. Works for both interior stencils (AbstractVector) and Hermite stencils (HermiteStencilData) via dispatch helpers in _build_collocation_matrix! and _build_rhs!.

Returns: view of the first k rows of λ (size k×num_ops).

source

RadialBasisFunctions._build_weights Method

_build_weights(ℒ, data, eval_points, adjl, basis)

Apply operator to basis functions and route to weight computation.

source

RadialBasisFunctions._build_weights Method

_build_weights(data, eval_points, adjl, basis, ℒrbf, ℒmon, mon;
              batch_size=10, device=CPU())

Build weights for interior-only problems (no boundary conditions). Creates empty BoundaryData to indicate all interior points.

source

RadialBasisFunctions._build_weights Method

_build_weights(ℒ, data, eval_points, adjl, basis,
              is_boundary, boundary_conditions, normals)

Generic Hermite dispatcher for operators. Applies operator to basis and routes to Hermite weight computation.

This eliminates repetitive _build_weights methods across operator files. Note: Type constraint removed to avoid circular dependency with operators.jl

source

RadialBasisFunctions._build_weights Method

_build_weights(data, eval_points, adjl, basis, ℒrbf, ℒmon, mon,
              is_boundary, boundary_conditions, normals;
              batch_size=10, device=CPU())

Build weights for problems with boundary conditions using Hermite interpolation. Exact allocation: Dirichlet points get single entry, others get full stencil.

source

RadialBasisFunctions._build_weights Method

_build_weights(ℒ, op)

Entry point from operator construction. Extracts configuration from operator and routes to appropriate implementation.

source

RadialBasisFunctions._expand_div_grad Method

_expand_div_grad(κ::AbstractVector)

Expand ∇⋅(κ∇) with per-dimension diffusivity into ∑ κ[i] * ∂²(i).

source

RadialBasisFunctions._expand_div_grad Method

_expand_div_grad(κ::Number)

Expand ∇⋅(κ∇) with scalar diffusivity into κ * Laplacian().

source

RadialBasisFunctions._expand_dot_grad Method

_expand_dot_grad(c::AbstractVector)

Expand c ⋅ ∇ into ∑ c[i] * ∂/∂xᵢ (advection operator).

source

RadialBasisFunctions._forward_with_cache Method

_forward_with_cache(data, eval_points, adjl, basis, ℒrbf, ℒmon, mon, ℒType)

Forward pass that builds weights while caching intermediate results for backward pass.

Returns: (W, cache) where W is the StencilWeights weight matrix and cache contains per-stencil factorizations and solutions needed for the pullback. Must produce bitwise the same weights as the primal build so the two paths are interchangeable.

source

RadialBasisFunctions._get_grad_funcs Method

_get_grad_funcs(OpType, basis, ℒ)

Get gradient functions for the given operator type and basis. Returns (grad_Lφ_x, grad_Lφ_xi) tuple.

source

RadialBasisFunctions._get_rhs_closures Method

_get_rhs_closures(OpType, ℒ, basis)

Get operator-specific closures for backward_stencil_with_ε!. Returns (poly_backward!, ∂Lφ_∂ε_fn) keyword arguments.

• Partial: polynomial section backward + partial ε derivative

• Laplacian: no polynomial backward + laplacian ε derivative

source

RadialBasisFunctions._hessian_row Method

_hessian_row(basis, x, xᵢ, dim)

Compute ∂/∂x[j] of [∂φ/∂x[dim]]: row dim of the Hessian of φ(x, xᵢ) w.r.t. x, extracted from the basis Hessian functor H(basis). x is a StaticVector (the only point type the backward pass produces); returns an SVector.

source

RadialBasisFunctions._hessian_row Method

_hessian_row(::PHS1, x, xᵢ, dim)

PHS1 needs exact handling: at r = 0 the derivatives of φ vanish so the gradient contribution is 0, but the raw H{PHS1} functor regularizes with avoid_inf (≈1e16 diagonals at r = 0, and an inexact r³ + 1e-16 denominator that distorts rows for r ≲ 1e-5), so the row is computed directly.

source

RadialBasisFunctions._infer_rank Method

_infer_rank(ℒ)

Infer the tensor rank of a custom operator function by probing it with a default basis. Returns 1 if ℒ(basis) produces a Tuple (one callable per dimension), 0 otherwise.

source

RadialBasisFunctions._interpolator_constructor_backward Method

_interpolator_constructor_backward(Δrbf_weights, Δmon_weights, A, k)

Backward pass for the Interpolator constructor w.r.t. y (the data values).

Given cotangents of rbf_weights and monomial_weights, computes the cotangent of y using the implicit function theorem. Since w = A⁻¹ [y; 0] and A is constant w.r.t. y:

Δy = (A⁻¹ [Δrbf_weights; Δmon_weights])[1:k]

A may be the collocation matrix or a factorization of it (both extension rules pass the cached BunchKaufman; symmetric ⟹ self-adjoint, so A⁻ᵀ = A⁻¹ holds either way).

Used by the Enzyme extension.

source

RadialBasisFunctions._interpolator_point_gradient! Method

_interpolator_point_gradient!(Δx, interp::Interpolator, x, Δy)

Accumulate the gradient of interp(x) * Δy into Δx.

RBF contribution: Σᵢ wᵢ * Δy * ∇φ(x, xᵢ) Polynomial contribution: Σⱼ wⱼ * Δy * ∇pⱼ(x)

source

RadialBasisFunctions._interpolator_weight_cotangents! Method

_interpolator_weight_cotangents!(Δrbf_w, Δmon_w, interp::Interpolator, x, Δy)

Accumulate the cotangent of interp(x) * Δy w.r.t. the interpolator weights:

Δrbf_w[i] += Δy * φ(x, xᵢ)
Δmon_w[j] += Δy * pⱼ(x)

Used by the Enzyme extension's Duplicated-Interpolator evaluation rules to deposit weight cotangents into the shadow Interpolator consumed by the constructor rule.

source

RadialBasisFunctions._optype Method

_optype(ℒ)

Map operator instance to its abstract type for dispatch in AD rules.

source

RadialBasisFunctions.accumulate_eval_pullback! Method

accumulate_eval_pullback!(Δx, W, Δy)

Accumulate the input cotangent of an operator application y = W * x into Δx: Δx += Wᵀ Δy. Used by the Enzyme eval rules. For StencilWeights this is the deterministic transpose-map adjoint gather (threaded on CPU, kernel on GPU; Dirichlet pad slots contribute zero); the generic method covers the sparse fallback (VirtualPartial). Matrix cotangents (operators applied to multi-column fields) accumulate column-wise.

source

RadialBasisFunctions.accumulate_weight_cotangent! Method

accumulate_weight_cotangent!(ΔWvals, W, x, Δy)

Accumulate the weight cotangent of y = W * x into the stencil-major value-matrix cotangent ΔWvals (the tangent of parent(W)): ΔW[l, i] += Δy[i] * x[idx[l, i]] — a pure gather over the stencil graph. Used by the Enzyme rule for * on StencilWeights when the weights themselves are active.

source

RadialBasisFunctions.allocate_ell Method

allocate_ell(backend, TD, k, N_eval, num_ops)

Allocate ELL (stencil-wise) weight storage on backend: one dense k × N_eval value matrix per operator component and a shared k × N_eval Int32 neighbor-index matrix.

source

RadialBasisFunctions.backward_collocation! Method

backward_collocation!(Δdata, ΔA, ∇p, neighbors, data, basis, mon, k)

Chain rule through collocation matrix construction.

The collocation matrix has structure: A[i,j] = φ(xi, xj) for i,j ≤ k (RBF block) A[i,k+j] = pⱼ(xi) for i ≤ k (polynomial block)

For RBF block (using ∇φ from existing basis_rules): Δxi += ΔA[i,j] * ∇φ(xi, xj) Δxj -= ΔA[i,j] * ∇φ(xi, xj) (by symmetry of φ(x-y))

For polynomial block: Δxi += ΔA[i,k+j] * ∇pⱼ(xi)

Note: A is symmetric, so we need to handle both triangles.

source

RadialBasisFunctions.backward_collocation_ε! Method

backward_collocation_ε!(Δε_acc, ΔA, neighbors, data, basis, k)

Compute gradient contribution to shape parameter ε from collocation matrix.

Uses implicit differentiation: Δε += Σᵢⱼ ΔA[i,j] * ∂A[i,j]/∂ε where A[i,j] = φ(xi, xj) for the RBF block.

source

RadialBasisFunctions.backward_linear_solve! Method

backward_linear_solve!(ΔA, Δb, Δλ, Δw, cache)

Compute cotangents of collocation matrix A and RHS vector b from cotangent of weights Δw. Δλ is a caller-owned scratch buffer (reused across stencils).

Given: Aλ = b, w = λ[1:k] We have: Δλ = [Δw; 0] (padded with zeros for monomial part)

Using implicit function theorem: η = A⁻ᵀ Δλ ΔA = -η λᵀ Δb = η

source

RadialBasisFunctions.backward_rhs! Method

backward_rhs!(Δdata, Δeval_point, Δb, neighbors, eval_point, data, basis, k, grad_Lφ_x, grad_Lφ_xi; poly_backward!=nothing)

Chain rule through RHS vector construction for any operator.

RBF section (shared by all operators): b[i] = ℒφ(eval_point, xi) → accumulate ∂/∂eval_point and ∂/∂xi

Polynomial section (Partial only — Laplacian gives constants, no gradient): b[k+j] = ℒpⱼ(eval_point) → passed as optional poly_backward! closure

source

RadialBasisFunctions.backward_rhs_ε! Method

backward_rhs_ε!(Δε_acc, Δb, neighbors, eval_point, data, basis, k, ∂Lφ_∂ε_fn)

Compute gradient contribution to shape parameter ε from RHS.

∂Lφ_∂ε_fn(x, xi) returns ∂(ℒφ)/∂ε for the specific operator:

• Laplacian: (x, xi) -> ∂Laplacian_φ_∂ε(basis, x, xi)

• Partial: (x, xi) -> ∂Partial_φ_∂ε(basis, dim, x, xi)

source

RadialBasisFunctions.backward_stencil_with_ε! Method

backward_stencil_with_ε!(Δdata, Δeval_point, Δε_acc, Δw, ΔA, Δb, Δλ, ∇p, cache, neighbors, eval_point, data, basis, mon, k, grad_Lφ_x, grad_Lφ_xi; poly_backward!=nothing, ∂Lφ_∂ε_fn=nothing)

Complete backward pass for a single stencil including shape parameter gradient.

Combines: 2. backward_linear_solve! → compute ΔA, Δb from Δw

3. backward_collocation! → chain ΔA to Δdata

4. backward_collocation_ε! → chain ΔA to Δε

5. backward_rhs! → chain Δb to Δdata and Δeval_point

6. backward_rhs_ε! → chain Δb to Δε

Optional closures:

• poly_backward!: polynomial section gradient (Partial only)

• ∂Lφ_∂ε_fn: shape parameter derivative function (IMQ/Gaussian only)

source

RadialBasisFunctions.build_weights_kernel Method

build_weights_kernel(data, eval_points, adjl, basis, ℒrbf, ℒmon, mon,
                    boundary_data; batch_size, device)

Main orchestrator for weight computation. Currently CPU-only. GPU stencil solve is not yet supported — see GitHub issue #88.

source

RadialBasisFunctions.build_weights_pullback_loop! Method

build_weights_pullback_loop!(Δdata, Δeval, Δε_acc, ΔW_extractor, ws, cache, adjl,
    eval_points, data, basis, mon, ℒ, OpType, grad_Lφ_x, grad_Lφ_xi)

Shared stencil iteration loop for _build_weights pullback across all AD backends.

ΔW_extractor(Δw, eval_idx, neighbors, k) is a callable that fills the stencil cotangent buffer Δw given the eval index, neighbor list, and stencil size. Both backends store cotangents as the stencil-major (k × N_eval) tangent of parent(::StencilWeights), so the extractors are column copies; the callable indirection remains so a backend with a different layout can still plug in.

The whole pullback path is host-pinned and stencil-major by design: the AD forward cache allocates plain Matrix{TD} weight buffers (mirroring the CPU-only primal build), this loop scalar-indexes them per stencil, and extract_stencil_cotangent! requires a host Matrix so a device array or device-oriented (SELL-resliced) layout cannot leak in silently. Device-adapted operators do not differentiate through their device weights — weight construction, and therefore its pullback, happens on host.

source

RadialBasisFunctions.calculate_batch_range Method

Calculate batch index range for kernel execution

source

RadialBasisFunctions.compute_hermite_poly_entry! Method

compute_hermite_poly_entry!(a, i, data, mon)

Compute polynomial entries for Hermite interpolation. Dispatches based on boundary point type for type stability.

source

RadialBasisFunctions.compute_hermite_rbf_entry Method

compute_hermite_rbf_entry(i, j, data, ops)

Compute single RBF matrix entry with Hermite boundary modifications. Dispatches based on point types (Interior/Dirichlet/NeumannRobin).

Uses BasisOperators for efficient evaluation (avoids functor construction in hot loop).

source

RadialBasisFunctions.construct_global_to_boundary Method

construct_global_to_boundary(is_boundary)

Construct mapping from global indices to boundary-only indices. For boundary points: global_to_boundary[i] = boundary array index For interior points: global_to_boundary[i] = 0 (sentinel)

source

RadialBasisFunctions.extract_stencil_cotangent! Method

extract_stencil_cotangent!(Δw, ΔW, eval_idx, k)

Fill the caller-owned buffer Δw with a single stencil's cotangent — column eval_idx of the stencil-major (k × N_eval) cotangent matrix ΔW (the tangent of parent(::StencilWeights)). Used by the Enzyme extension.

ΔW is deliberately constrained to a host Matrix: the AD path is structurally host-pinned (the forward cache allocates Matrix{TD} and the pullback loop scalar-indexes), so a device array or a device-oriented (resliced) values layout leaking in here fails with a MethodError instead of silently mis-indexing.

source

RadialBasisFunctions.fill_dirichlet_column! Method

Fill a Dirichlet identity column: weight 1 at slot 1, zero pads elsewhere. All slots carry the eval point's own index (in-bounds; Dirichlet stencils require eval_points === data), so duplicate-combining sparse conversion collapses the column to a single identity entry.

source

RadialBasisFunctions.fill_entries! Method

Write one stencil's weights into its ELL column

source

RadialBasisFunctions.grad_applied_laplacian_wrt_x Method

grad_applied_laplacian_wrt_x(basis)

Get gradient of applied Laplacian operator w.r.t. evaluation point.

source

RadialBasisFunctions.grad_applied_laplacian_wrt_xi Method

grad_applied_laplacian_wrt_xi(basis)

Get gradient of applied Laplacian operator w.r.t. data point. By symmetry, always the negation of the _wrt_x version.

source

RadialBasisFunctions.grad_applied_partial_wrt_x Method

grad_applied_partial_wrt_x(basis, dim)

Get gradient of applied partial derivative operator w.r.t. evaluation point.

source

RadialBasisFunctions.grad_applied_partial_wrt_xi Method

grad_applied_partial_wrt_xi(basis, dim)

Get gradient of applied partial derivative operator w.r.t. data point. By symmetry, always the negation of the _wrt_x version.

source

RadialBasisFunctions.grad_laplacian_gaussian_wrt_x Method

grad_laplacian_gaussian_wrt_x(ε)

Returns a function computing ∂/∂x[j] of [∇²φ] for Gaussian.

Mathematical derivation: ∇²φ = (4ε⁴r² - 2ε²D) * φ

∂(∇²φ)/∂x[j] = φ * δ_j * 4ε⁴ * [2 + D - 2ε²r²]

source

RadialBasisFunctions.grad_laplacian_imq_wrt_x Method

grad_laplacian_imq_wrt_x(ε)

Returns a function computing ∂/∂x[j] of [∇²φ] for IMQ.

Mathematical derivation: ∇²φ = sum_i [∂²φ/∂x[i]²] = 3ε⁴r²/s^(5/2) - D*ε²/s^(3/2)

∂(∇²φ)/∂x[j] = δ_j * [3(D+2)ε⁴/s^(5/2) - 15ε⁶r²/s^(7/2)]

source

RadialBasisFunctions.grad_laplacian_phs1_wrt_x Method

grad_laplacian_phs1_wrt_x()

Returns a function computing ∂/∂x[j] of [∇²φ] for PHS1.

Mathematical derivation: ∇²φ = (d−1)/r ∂/∂x[j] [(d−1)/r] = -(d−1) * δ_j / r³

Note: At r=0, we return 0 to avoid singularity.

source

RadialBasisFunctions.grad_laplacian_phs3_wrt_x Method

grad_laplacian_phs3_wrt_x()

Returns a function computing ∂/∂x[j] of [∇²φ] for PHS3.

Mathematical derivation: ∇²φ = 3(d+1)r ∂/∂x[j] [3(d+1)r] = 3(d+1) * δ_j / r

source

RadialBasisFunctions.grad_laplacian_phs5_wrt_x Method

grad_laplacian_phs5_wrt_x()

Returns a function computing ∂/∂x[j] of [∇²φ] for PHS5.

Mathematical derivation: ∇²φ = 5(d+3)r³ ∂/∂x[j] [5(d+3)r³] = 5(d+3) * 3r² * δ_j / r = 15(d+3) * r * δ_j

source

RadialBasisFunctions.grad_laplacian_phs7_wrt_x Method

grad_laplacian_phs7_wrt_x()

Returns a function computing ∂/∂x[j] of [∇²φ] for PHS7.

Mathematical derivation: ∇²φ = 7(d+5)r⁵ ∂/∂x[j] [7(d+5)r⁵] = 7(d+5) * 5r⁴ * δ_j / r = 35(d+5) * r³ * δ_j

source

RadialBasisFunctions.hermite_mono_rhs! Method

hermite_mono_rhs!(bmono, ℒmon, mon, eval_point, is_bound, bc, normal,
                  workspace, normal_workspace)

Apply boundary conditions to monomial operator evaluation. Dispatches based on boundary point type for type stability.

source

RadialBasisFunctions.hermite_rbf_rhs Method

hermite_rbf_rhs(ℒrbf, eval_point, data_point, is_bound, bc, normal)

Apply boundary conditions to RBF operator evaluation.

• Interior/Dirichlet: standard evaluation ℒΦ(x_eval, x_data)

• Neumann/Robin: α_ℒΦ + β_ℒ(∂ₙΦ)

source

RadialBasisFunctions.launch_kernel! Method

launch_kernel!(...)

Launch parallel CPU kernel for weight computation. Handles Dirichlet/Interior/Hermite stencil classification via dispatch.

source

RadialBasisFunctions.negate_grad Method

negate_grad(grad_fn)

Given a gradient function grad_fn(x, xi), returns (x, xi) -> -grad_fn(x, xi). All _wrt_xi functions are the negation of their _wrt_x counterparts by symmetry.

source

RadialBasisFunctions.point_type Method

Determine boundary type of a single point

source

RadialBasisFunctions.run_build_weights_pullback Method

run_build_weights_pullback(ΔW_extractor, cache, adjl, eval_points, data, basis,
    mon, ℒ, OpType, grad_Lφ_x, grad_Lφ_xi)

Allocate cotangent buffers, run build_weights_pullback_loop!, and return (Δdata, Δeval, Δε_acc). Used by the Enzyme extension.

source

RadialBasisFunctions.∂Laplacian_φ_∂ε Method

∂Laplacian_φ_∂ε(basis::Gaussian, x, xi)

Derivative of Laplacian of Gaussian basis w.r.t. shape parameter ε.

∇²φ = (-2ε²D + 4ε⁴r²) exp(-ε²r²), where D = dimension ∂(∇²φ)/∂ε = exp(-ε²r²) [-4εD + 16ε³r² + 4ε³r²D - 8ε⁵r⁴]

source

RadialBasisFunctions.∂Laplacian_φ_∂ε Method

∂Laplacian_φ_∂ε(basis::IMQ, x, xi)

Derivative of Laplacian of IMQ basis w.r.t. shape parameter ε.

Let s = ε²r² + 1, D = dimension ∇²φ = -ε²D/s^(3/2) + 3ε⁴r²/s^(5/2) ∂(∇²φ)/∂ε = ∂/∂ε[-ε²D s^(-3/2) + 3ε⁴r² s^(-5/2)]

source

RadialBasisFunctions.∂Partial_φ_∂ε Method

∂Partial_φ_∂ε(basis::Gaussian, dim::Int, x, xi)

Derivative of first partial derivative of Gaussian basis w.r.t. shape parameter ε.

∂φ/∂x_dim = -2ε²(x_dim - xi_dim) exp(-ε²r²) ∂/∂ε[∂φ/∂x_dim] = 4ε(x_dim - xi_dim)(ε²r² - 1) exp(-ε²r²)

source

RadialBasisFunctions.∂Partial_φ_∂ε Method

∂Partial_φ_∂ε(basis::IMQ, dim::Int, x, xi)

Derivative of first partial derivative of IMQ basis w.r.t. shape parameter ε.

∂φ/∂x_dim = ε²(xi_dim - x_dim) s^(-3/2) ∂/∂ε[∂φ/∂x_dim] = 2ε(xi_dim - x_dim) s^(-3/2) + ε²(xi_dim - x_dim)(-3/2)s^(-5/2) · 2εr² = (xi_dim - x_dim)[2ε s^(-3/2) - 3ε³r² s^(-5/2)]

source

RadialBasisFunctions.∂_normal! Method

∂_normal!(b, scratch, mb, normal, x)

In-place normal derivative of the monomial basis: b = Σ_d normal[d] * ∂_d P(x). scratch must have the same length as b.

source

RadialBasisFunctions.∂φ_∂ε Method

∂φ_∂ε(basis::Gaussian, x, xi)

Derivative of Gaussian basis function w.r.t. shape parameter ε.

φ(r) = exp(-ε²r²) ∂φ/∂ε = -2εr² exp(-ε²r²)

source

RadialBasisFunctions.∂φ_∂ε Method

∂φ_∂ε(basis::IMQ, x, xi)

Derivative of IMQ basis function w.r.t. shape parameter ε.

φ(r) = (ε²r² + 1)^(-1/2) ∂φ/∂ε = -εr² (ε²r² + 1)^(-3/2)

source

SparseArrays.sparse Method

sparse(op::RadialBasisOperator)
SparseMatrixCSC(op::RadialBasisOperator)

Convert the operator's weights to a SparseMatrixCSC — or a tuple of them for gradient-family operators — for global PDE system assembly or interop with sparse linear algebra libraries. Rebuilds the weights first if the cache is stale.

source

RadialBasisFunctions.AbstractJacobianOperator Type

AbstractJacobianOperator{Dim,N} <: AbstractOperator{N}

Abstract supertype for operators whose action on a basis is the tuple of first partial derivatives (∂₁, …, ∂_Dim): Jacobian, Divergence, Curl, StrainRate, and RotationRate.

source

RadialBasisFunctions.BackwardWorkspace Type

BackwardWorkspace{T,P,M,V}

Per-pass scratch for the _build_weights backward pass, allocated once and reused across every stencil to eliminate per-stencil heap allocations. All buffer shapes are fixed and uniform for a pass (n = k + nmon, num_ops, dim_space), so they are sized once from the forward cache.

Container fields are type-parameterized (M, V) rather than hard-coded Matrix/Vector, matching StencilForwardCache; the constructor builds concrete CPU arrays and lets the parameters be inferred.

• ΔA: n × n adjoint of the collocation matrix

• Δb: n × num_ops adjoint of the RHS

• Δλ: n × num_ops padded stencil cotangent (overwritten by the adjoint solve)

• ∇p: nmon × dim_space monomial-gradient scratch

• Δw: k × num_ops extracted stencil cotangent

• Δlocal_data: k buffers of length dim_space (local data-point cotangents)

• Δeval_pt: length dim_space eval-point cotangent

• local_data: k gathered stencil points

• local_idx: [1, …, k], the constant local neighbor indices

source

RadialBasisFunctions.BasisOperators Type

BasisOperators{B,G,Hess}

Bundle of pre-constructed basis operators for efficient evaluation in hot loops. Avoids repeated functor construction inside hermite_rbf_dispatch.

Fields:

• φ: The basis function itself

• ∇φ: Gradient operator (pre-constructed ∇(basis))

• Hφ: Hessian operator (pre-constructed H(basis))

Usage:

ops = BasisOperators(basis)
# In hot loop:
φ_val = ops.φ(x, xᵢ)
grad = ops.∇φ(x, xᵢ)      # Returns vector
hess = ops.Hφ(x, xᵢ)      # Returns matrix
Dφ = dot(n, grad)         # Directional derivative
D²φ = dot(ni, hess * nj)  # Second directional derivative

source

RadialBasisFunctions.BasisOperators Method

Construct BasisOperators from a basis function.

source

RadialBasisFunctions.BoundaryData Type

BoundaryData{T,V}

Wrapper for global boundary information (replaces fragile tuples).

source

RadialBasisFunctions.BoundaryPointType Type

Trait types for individual point boundary classification

source

RadialBasisFunctions.H Type

H{B<:AbstractRadialBasis}

Hessian operator functor. Construct with H(basis). Returns the Hessian matrix of the basis function.

source

RadialBasisFunctions.ScaledKernel Type

ScaledKernel{T<:Number, F}

Callable scaling a basis action by α. Supports the standard (x, xᵢ) form and the Hermite normal form (x, xᵢ, normal).

source

RadialBasisFunctions.StencilForwardCache Type

StencilForwardCache{T}

Per-stencil storage from forward pass needed for backward pass.

• lambda: Full solution vector (k+nmon) × num_ops from solving Aλ = b

• A_fact: Bunch-Kaufman factorization of the symmetric collocation matrix

• k: Number of RBF neighbors in stencil

• nmon: Number of monomial basis functions

source

RadialBasisFunctions.StencilType Type

Trait types for stencil classification

source

RadialBasisFunctions.SumKernel Type

SumKernel{Fs<:Tuple}

Callable summing component basis actions. Supports the standard (x, xᵢ) form and the Hermite normal form (x, xᵢ, normal).

source

RadialBasisFunctions.WeightsBuildForwardCache Type

WeightsBuildForwardCache{T}

Global cache storing all stencil results and references to inputs.

• stencil_caches: Vector of StencilForwardCache, one per evaluation point

• k: Stencil size (number of neighbors)

• nmon: Number of monomial basis functions

• num_ops: Number of operators (1 for scalar, D for gradient)

source

RadialBasisFunctions.ℒMonomialBasis Type

ℒMonomialBasis{Dim,Deg,F}

A differentiated MonomialBasis: wraps an in-place evaluator f(b, x) that fills b with the operator applied to each monomial term at x.

This is the monomial half of the package's two basis-differentiation protocols: for MonomialBasis, ∂/∇/∂²/∇²/H/∂mixed are factory functions returning an ℒMonomialBasis, whereas for AbstractRadialBasis the same names are functor structs (see src/basis/basis.jl). Operator actions such as (op::Partial)(basis) are written once and serve both protocols.

source

RadialBasisFunctions.∂ Type

∂{B<:AbstractRadialBasis}

Partial derivative operator functor. Construct with ∂(basis, dim).

source

RadialBasisFunctions.∂mixed Type

∂mixed{B<:AbstractRadialBasis}

Mixed partial derivative operator functor. Construct with ∂mixed(basis, dim1, dim2). Computes ∂²φ/(∂x_{dim1} ∂x_{dim2}).

source

RadialBasisFunctions.∂² Type

∂²{B<:AbstractRadialBasis}

Second partial derivative operator functor. Construct with ∂²(basis, dim).

source

RadialBasisFunctions.∇ Type

∇{B<:AbstractRadialBasis}

Gradient operator functor. Construct with ∇(basis).

source

RadialBasisFunctions.∇² Type

∇²{B<:AbstractRadialBasis}

Laplacian operator functor. Construct with ∇²(basis).

source

RadialBasisFunctions.VectorFieldOperator Type

requires_vector_input(op) -> Bool

Whether the operator requires a vector field (matrix) as input rather than a scalar field (vector). Operators like Divergence, Curl, StrainRate, and RotationRate act on vector fields and will error on scalar input.

source

EllSparse

Self-contained SELL-C/ELL sparse storage backing the stencil weight matrices. Not re-exported — reach these names as RadialBasisFunctions.EllSparse.<name>. See the internals for how it fits the weight pipeline.

RadialBasisFunctions.EllSparse Module

EllSparse

Self-contained SELL-C (Sliced ELLpack) sparse matrix storage with CPU and GPU KernelAbstractions kernels. SellMatrix{T, C} stores an m × n sparse matrix as horizontal slices of C consecutive rows, each slice padded to its own width:

         slots (l) →
row 1  ┌ a  b  c ┐
row 2  └ d  e  ⋅ ┘   slice 1 (C = 2 rows, width w₁ = 3; ⋅ = sentinel padding)
row 3  ┌ f  ⋅    ┐
row 4  └ g  h    ┘   slice 2 (width w₂ = 2)

Within slice s the entry for (local row r, slot l) lives at flat storage position

q = sliceptr[s] + (l - 1) * C + (r - 1)

i.e. slot-major within a slice: for a fixed slot, consecutive rows are adjacent in memory. With C matched to the device warp/wavefront size (32), a row-parallel kernel reads coalesced — this is the cuSPARSE SELL layout (modulo their -1 padding sentinel vs 0 here, and 0- vs 1-based offsets). With C = 1 the layout degenerates to row-contiguous ELL: q = rowstart + l - 1, and for a uniform width k the flat storage is bit-identical to a column-major k × m matrix whose column i holds row i's entries — the natural cache-friendly CPU layout for stencil-style matrices.

Format decisions

• One type covers both orientations: slice height C lives in the type domain (SellMatrix{T, C, ...}), with EllMatrix = SellMatrix{T, 1} as an alias. Layout changes are explicit via reslice; Adapt.adapt never changes layout (a data move must not change adjoint summation order).

• Classic single-slice slot-major ELL is intentionally absent — it was only ever a GPU format (slot-major striding is cache-hostile for CPU row applies, where C = 1 wins), and on GPUs C = 32 strictly dominates it (cuSPARSE itself dropped plain ELL in favor of SELL). Its global-max-width padding survives as the constructor option pad = :global.

• Structure and values are split: everything except vals lives in an immutable SellStructure that multiple matrices alias (===) via with_values, so families of same-pattern matrices (e.g. gradient components) share one column-index array and one transpose map.

• Adjoint/transpose products are deterministic and atomics-free: a precomputed TransposeMap turns them into per-column gathers with a fixed summation order (ascending storage position), stable under any thread count.

Scoped-out extension points

• SELL-C-σ row sorting (a row permutation before slicing; zero benefit at uniform row length and it complicates AD seams).

• COO ingestion (construct via SparseMatrixCSC for now).

• Asynchronous batched mul! (multi-column applies enqueue all columns, then synchronize once).

• cuSPARSE interop: their SELL uses sentinel -1 and 0-based offsets; a future CUDA extension converts at the boundary.

• Parallel transpose-map build (a chunked counting sort with per-chunk cursors can preserve the ascending-position order, but the build is a once-per-structure host pass dwarfed by weight computation — not worth risking the order contract).

source

RadialBasisFunctions.EllSparse.EllMatrix Type

EllMatrix{T, V, S}

Alias for SellMatrix{T, 1, V, S}: slice height 1, i.e. row-contiguous ELL storage. For a uniform stored width k the flat storage is bit-identical to a column-major k × m matrix of per-row entries — the cache-friendly CPU orientation.

source

RadialBasisFunctions.EllSparse.SellMatrix Type

SellMatrix{T, C, V, S} <: AbstractMatrix{T}

SELL-C sparse matrix: a values array over an immutable SellStructure. See the EllSparse module docstring for the storage layout. A * x computes y[i] = Σₗ vals[q(i, l)] * x[colind[q(i, l)]] row-parallel on CPU threads or GPU; A' * x is a deterministic per-column gather through the structure's transpose map.

vals is a flat vector of length sliceptr[end] - 1, except in the uniform C == 1 case where it may be the k × m matrix the layout is bit-identical to (zero-copy: parent(A) === vals).

source

RadialBasisFunctions.EllSparse.SellMatrix Method

SellMatrix(vals::AbstractMatrix, colind::AbstractMatrix{<:Integer}, n::Integer;
           transpose_map = true)

Zero-copy uniform C == 1 construction from per-row matrices: vals[l, i] is the l-th stored value of row i and colind[l, i] its column, both k × m (column per logical row). The flat storage aliases both inputs (parent(A) === vals); every slot must hold a real entry (colind values in 1:n, no sentinels).

Set transpose_map = false to skip building the TransposeMap; adjoint and transpose products then throw until the matrix is rebuilt with a map.

source

RadialBasisFunctions.EllSparse.SellMatrix Method

SellMatrix(S::SparseMatrixCSC, Val(C); pad = :slice, index_type = Int32,
           transpose_map = true)

Convert a CSC matrix to SELL-C. Row lengths are counted, each slice is padded to the width of its longest row (pad = :slice) or to the global maximum (pad = :global), and entries are filled in CSC column order so each row's column indices come out ascending. Explicit stored zeros are retained. Host-serial; device conversion goes through a host matrix.

See also sell for a runtime slice height.

source

RadialBasisFunctions.EllSparse.SellStructure Type

SellStructure{C, Ti, VI, TM}

The immutable sparsity structure of a SellMatrix with slice height C: logical size, column indices, slice offsets, and the optional TransposeMap. Multiple matrices alias one structure (===) via with_values, which makes same-pattern families (shared column indices, shared transpose map) first-class.

colind is the flat column-index array; the sentinel Ti(0) marks padding slots, including the ghost rows of a last partial slice when m % C != 0. The values array of any matrix over this structure stores zero(T) in every padding slot — an invariant the same-structure ==/isapprox fast paths rely on. width is the uniform slice width when every slice has the same width, else 0. padded is true when any sentinel exists and gates the unguarded (branchless) kernels.

source

RadialBasisFunctions.EllSparse.TransposeMap Type

TransposeMap{Ti, VI, VR}

Column-wise adjacency of a SELL structure: for column j, positions[offsets[j]:offsets[j+1]-1] are the flat storage positions of the stored entries whose column index is j. The per-column sequence is the adjoint summation order — the determinism anchor: adjoint/transpose products walk this map as a per-column gather with a fixed order, so results are bitwise-reproducible under any thread count, with no atomics. build_transpose_map produces ascending positions; reslice remaps them preserving the sequence (not re-sorting), so a layout change never changes the summation order.

rows carries the global row of each entry (aligned with positions) and is nothing exactly when C == 1 with a uniform stored width k, where the row is recoverable as (q - 1) ÷ k + 1.

source

Base.parent Method

parent(A::SellMatrix)

The raw values storage backing A (zero-copy): the flat vector, or the k × m matrix in the uniform C == 1 matrix-backed case. The supported handle for reading or mutating stored values in place.

source

RadialBasisFunctions.EllSparse.adapt_family Method

adapt_family(to, family::NTuple{N, SellMatrix}) -> NTuple{N, SellMatrix}

Adapt a family of same-structure matrices (e.g. gradient components) to a backend, adapting the shared structure once and rebinding each member's values onto it — so the === structure aliasing survives the move. All members must alias one structure object.

source

RadialBasisFunctions.EllSparse.build_transpose_map Method

build_transpose_map(Val(C), colind, sliceptr, m, n, width) -> TransposeMap

Build the column-wise TransposeMap of a SELL structure by counting sort: one linear pass over the flat colind (skipping sentinels) counts entries per column, a cumulative sum forms offsets, and a second linear pass in ascending storage position fills positions — which is what makes adjoint summation order deterministic. rows is filled alongside unless C == 1 with uniform width > 0, where the row is recoverable as (q - 1) ÷ width + 1.

The sort is a serial host pass: device inputs are copied to the host once and the finished map is copied back, so the map lives uniformly on colind's backend. Never scalar-indexes device arrays.

source

RadialBasisFunctions.EllSparse.coo Method

coo(A::SellMatrix) -> (I, J, V)

The stored (non-sentinel) entries of A as COO triplets, in row order with each row's slots ascending. Explicit stored zeros are retained; duplicates are not combined here — sparse(A) sums them.

source

RadialBasisFunctions.EllSparse.preferred_slice_height Method

preferred_slice_height(backend) -> Val

The slice height a backend's row-parallel kernels want: Val(1) on CPU (row-contiguous SIMD), Val(32) on device backends (warp-coalesced, the cuSPARSE SELL orientation). The policy hook for layout decisions at data-movement boundaries; pair with reslice — Adapt.adapt itself never changes layout.

source

RadialBasisFunctions.EllSparse.reslice Method

reslice(A::SellMatrix, Val(C₂); pad = :slice) -> SellMatrix

Rebuild A with slice height C₂ — the explicit layout-change operation (adapt never changes layout). Stored entries keep their per-row slot order, and the transpose map is remapped in sequence order rather than rebuilt, so adjoint products of the resliced matrix keep the original summation order bit-for-bit (up to hardware FMA differences across devices).

Device matrices round-trip through the host (the map remap is host-serial regardless); the result lives on A's backend.

source

RadialBasisFunctions.EllSparse.same_structure Method

same_structure(A::SellMatrix, B::SellMatrix) -> Bool

Whether A and B share one sparsity structure: identical objects (===), or equal slice height, size, slice offsets, and column indices.

source

RadialBasisFunctions.EllSparse.sell Method

sell(S::SparseMatrixCSC; slice_height = 1, kwargs...) -> SellMatrix

Runtime-slice-height entry to the CSC constructor, through a Val function barrier so downstream code still specializes on C. Keyword arguments are forwarded to SellMatrix(S, Val(C); ...).

source

RadialBasisFunctions.EllSparse.slice_height Method

slice_height(A::SellMatrix) -> Int

The slice height C of the storage layout: 1 is the row-contiguous CPU orientation, 32 the coalesced device orientation. Change it with reslice.

source

RadialBasisFunctions.EllSparse.structure Method

structure(A::SellMatrix) -> SellStructure

The immutable sparsity structure backing A. Matrices from the same family alias one structure object (compare with === or same_structure).

source

RadialBasisFunctions.EllSparse.uniform_width Method

uniform_width(A::SellMatrix) -> Int

The uniform stored width when every slice has the same width, else 0.

source

RadialBasisFunctions.EllSparse.values_matrix Method

values_matrix(A::SellMatrix{T, 1}) -> AbstractMatrix{T}

The k × m per-row values matrix of a uniform-width C == 1 matrix — A.vals itself when matrix-backed, else a zero-copy reshape of the flat storage. Errors for C != 1 or ragged widths, where no such matrix view exists.

source

RadialBasisFunctions.EllSparse.with_values Method

with_values(A::SellMatrix, vals::AbstractVecOrMat) -> SellMatrix

A new matrix over the same structure object (===) with a different values array. vals must match the storage length and live on the same backend. Padding slots of vals must hold zero(eltype(vals)) — the same-structure == fast path relies on it.

source