From 8e654cfdc447f8b3a04461317ce310ec2fb438a1 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 02:46:56 +0100 Subject: [PATCH 01/12] Caller-held buffers and workspaces for the hot evaluation, solve and projection paths MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit evaluate!(buf, b, û, x, d), evaluate!(buf, s, x, d) and evaluate!(bufs, B, û, x, d) evaluate a spline in buffers that the caller holds. mass_workspace and mass_solve!(y, op, x, work) solve with a KroneckerMass; contract_workspace, contract! and l2_projection!(û, q, F, work) contract and project on a TensorProductQuadrature. Each reads 0 B, also on a basis whose axes differ in type, which the loops over the axes now unroll. The allocating forms build their buffers and call these, and return bit-for-bit the values of cb82e6f. Records K6: a periodic basis gives NaN derivatives to a ForwardDiff.Dual at both ends. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 38 ++++ KNOWN_ISSUES.md | 14 ++ docs/src/library.md | 4 + src/SimpleSplines.jl | 5 +- src/bspline.jl | 54 ++++- src/spline.jl | 4 + src/tensorproduct.jl | 313 +++++++++++++++++++++++------ test/Project.toml | 2 + test/integration/workspaces.jl | 347 +++++++++++++++++++++++++++++++++ test/quality/jet.jl | 25 +++ test/runtests.jl | 1 + 11 files changed, 740 insertions(+), 67 deletions(-) create mode 100644 test/integration/workspaces.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index e6b13fa..9f4f8ea 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,44 @@ All notable changes to SimpleSplines.jl are documented here. The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and the project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] — targeting 0.3.5 + +### New Features + +The hot evaluation, solve and projection paths take buffers that the caller holds, and then +allocate nothing. `evaluate!(buf, b, û, x, d)` and `evaluate!(buf, s, x, d)` evaluate a spline +on a one-dimensional basis with a buffer of `local_width(b)` entries; `evaluate!(bufs, B, û, x, +d)` does it on a `TensorProductBasis` with one buffer per axis. `mass_workspace(op, T)` and +`mass_solve!(y, op, x, work)` solve with a `KroneckerMass`; `contract_workspace(q, T)`, +`contract!(L, q, F, work, d)` and `l2_projection!(û, q, F, work)` contract and project on a +`TensorProductQuadrature`. Each buffer and workspace belongs to one task. The allocating forms +stay, and build their buffers and call these. + +Bytes per call, measured with `@allocated` after a warm call through a function barrier, in a +cold process on Julia 1.13.3, on `cb82e6f` and on this release. The one-dimensional basis is +clamped and cubic on 8 cells; the two-dimensional one is periodic × clamped, cubic, on 8 × 8 +cells. Each pair is `Float64` / `Float32`. + +| call | `cb82e6f` | the allocating form now | the caller-buffer form | +|:--|--:|--:|--:| +| `evaluate(b, û, x, d)` | 96 / 80 | 96 / 80 | 0 / 0 | +| `evaluate(B, û, x, d)` | 192 / 160 | 192 / 160 | 0 / 0 | +| `mass_solve!(y, op, x)` | 1232 / 1104 | 544 / 416 | 0 / 0 | +| `contract(q, F, d)` | 24064 / 13536 | 30128 / 15328 | 0 / 0 | +| `l2_projection!(û, q, F)` | 25296 / 14640 | 29344 / 14880 | 0 / 0 | + +The solve along the axes of a `KroneckerMass` and the weighting and contraction along the axes +of `contract` are unrolled over the tuple of axes, so an operator or a quadrature whose axes have +different types no longer dispatches dynamically: that is the drop of the three-argument +`mass_solve!` on the periodic × clamped basis. The allocating `contract` and `l2_projection!` +allocate more than before, because they build the whole workspace: two buffers of the largest +intermediate array, and the buffers of the mass solve. `mass_solve!` on a `KroneckerMass` checks +the length of `y` before it copies `x`, and raises a `DimensionMismatch` for a `y` of the wrong +length. The results of the allocating forms are bit-for-bit those of `cb82e6f`. + +`ForwardDiff` joins the test dependencies: the tests evaluate with `ForwardDiff.Dual` arguments +and buffers, which also allocate nothing. + ## [0.3.4] — 2026-10-31 ### Bug Fixes diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 424c4ac..3fda417 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -71,3 +71,17 @@ `Float32` spans throw on #39 and no `Float64` span does. - **kind:** defect - **found:** 2026-10-31; older than #39 + +### K6 · A periodic basis gives `NaN` derivatives to a `ForwardDiff.Dual` at `a` and at `b` + +- **location:** `src/bspline.jl`, `_reduce_argument(b::PeriodicBSplineBasis, x)` +- **evidence:** the argument is reduced as `a + mod(x - a, L)`, and ForwardDiff's `mod` gives a + `NaN` partial where `x - a` is a multiple of `L`: `mod(Dual{Nothing}(1.0, 1.0), 1.0)` is + `Dual{Nothing}(0.0,NaN)`. With `b = BSplineBasis(UniformMesh(8, 0.0 .. 1.0), 3, Periodic())` + and `û = collect(1.0:8)` on `cb82e6f`, `evaluate(b, û, Dual{Nothing}(1.0, 1.0))` and + `evaluate(b, û, Dual{Nothing}(0.0, 1.0))` both give `Dual{Nothing}(6.999999999999999,NaN)`, + while `Dual{Nothing}(prevfloat(1.0), 1.0)` gives `Dual{Nothing}(7.0,8.000000000000014)`. + `evaluate_all(b, Dual{Nothing}(1.0, 1.0))` gives `NaN` partials in all four entries. The values + are right; only the derivatives are lost. +- **kind:** defect +- **found:** 2026-11-02 diff --git a/docs/src/library.md b/docs/src/library.md index 291b590..3bfb33b 100644 --- a/docs/src/library.md +++ b/docs/src/library.md @@ -86,6 +86,7 @@ polynomial_reproduction ```@docs evaluate +evaluate! evaluate_all evaluate_all! BSplineDerivative @@ -124,6 +125,7 @@ KroneckerMass mass_operator mass_factors mass_solve! +mass_workspace ``` ## Tensor products @@ -137,6 +139,8 @@ quadratures quadrature_grid_size quadrature_sample contract +contract! +contract_workspace ``` ## Polar splines diff --git a/src/SimpleSplines.jl b/src/SimpleSplines.jl index 46720de..64c3cf0 100644 --- a/src/SimpleSplines.jl +++ b/src/SimpleSplines.jl @@ -48,7 +48,7 @@ include("boundary.jl") export AbstractBSplineBasis, BSplineBasis, PeriodicBSplineBasis export BSplineDerivative, PeriodicBSplineDerivative -export knotvector, evaluate, evaluate_all, evaluate_all!, mesh +export knotvector, evaluate, evaluate!, evaluate_all, evaluate_all!, mesh export boundary, basis_index, findcell, local_width, polynomial_reproduction include("bspline.jl") @@ -69,7 +69,8 @@ export mass_matrix, mass_factorization, stiffness_matrix, derivative_matrix, include("quadrature.jl") export TensorProductBasis, TensorProductQuadrature, KroneckerMass, ⊗ -export bases, quadratures, mass_factors, contract +export bases, quadratures, mass_factors, mass_workspace, contract, contract!, + contract_workspace export quadrature_sample, quadrature_grid_size include("tensorproduct.jl") diff --git a/src/bspline.jl b/src/bspline.jl index 7562166..373226d 100644 --- a/src/bspline.jl +++ b/src/bspline.jl @@ -625,22 +625,55 @@ end function evaluate(b::AbstractBSplineBasis{T}, û::AbstractVector{S}, x::Number, d::Integer = 0) where {T, S} + R = _evaltype(promote_type(T, S), typeof(x)) + evaluate!(Vector{R}(undef, local_width(b)), b, û, x, d) +end + +@doc raw""" + evaluate!(buf, b::AbstractBSplineBasis, û::AbstractVector, x::Number, d = 0) + evaluate!(buf, s::Spline, x::Number, d = 0) + +The value of [`evaluate`](@ref)`(b, û, x, d)`, or of `evaluate(s, x, d)`, computed in the +caller's buffer `buf`, which allocates nothing. + +`buf` must have [`local_width`](@ref)`(b)` entries, or a `DimensionMismatch` is raised; its +contents on entry do not matter, and on return they are the local block of +[`evaluate_all!`](@ref). Its element type is the type the sum is accumulated in, so it must +hold `promote_type(eltype(b), eltype(û), typeof(x))`: a `ForwardDiff.Dual` argument takes a +buffer of `Dual`s. + +A buffer belongs to one task. Two tasks that evaluate at once need a buffer each. + +```jldoctest +julia> b = BSplineBasis(UniformMesh(8, 0 .. 1), 3); û = collect(1.0:nbasis(b)); + +julia> buf = zeros(local_width(b)); + +julia> evaluate!(buf, b, û, 0.3, 1) == evaluate(b, û, 0.3, 1) +true +``` +""" +function evaluate!(buf::AbstractVector, b::AbstractBSplineBasis, û::AbstractVector, + x::Number, d::Integer = 0) length(û) == nbasis(b) || throw(DimensionMismatch( "the coefficient vector has $(length(û)) entries but the basis has $(nbasis(b))")) - R = _evaltype(promote_type(T, S), typeof(x)) - _evaluate_block!(Vector{R}(undef, local_width(b)), b, û, x, d) + w = local_width(b) + length(buf) == w || throw(DimensionMismatch( + "evaluate! needs a buffer of $(w) entries for this basis, got $(length(buf)); " * + "see `local_width`")) + _evaluate_block!(buf, b, û, x, d) end # The local-block sum, against a caller-supplied buffer. Only the block contributes, so this # is O(p) per point where summing `evaluate` over `eachindex(û)` is O(N). # -# Split out so that a sweep over many points allocates one buffer for the sweep rather than -# one per point. The buffer is the whole of the difference between this and BSplineKit at -# small `n`, which reaches the block through the spline order carried in its *type* and can -# therefore keep it on the stack; here the degree is a field, so the length is a run-time -# value. Passing a `Val`-sized `MVector` instead was measured and does not help -- it halves -# the allocation and returns the saving in dispatch, because `evaluate_all!` takes an -# `AbstractVector` and the buffer escapes into it either way. +# The buffer is the whole of the difference between this and BSplineKit at small `n`, which +# reaches the block through the spline order carried in its *type* and can therefore keep it +# on the stack; here the degree is a field, so the length is a run-time value. Passing a +# `Val`-sized `MVector` instead was measured and does not help -- it halves the allocation and +# returns the saving in dispatch, because `evaluate_all!` takes an `AbstractVector` and the +# buffer escapes into it either way. The allocation-free path is therefore a buffer that the +# caller holds, through `evaluate!`; a sweep over many points below holds one for the sweep. function _evaluate_block!(values::AbstractVector{R}, b::AbstractBSplineBasis, û::AbstractVector, x::Number, d::Integer) where {R} j₀ = evaluate_all!(values, b, x, d) @@ -883,7 +916,8 @@ The allocating form of [`evaluate_all!`](@ref): returns `(j₀, values)` with a buffer of [`local_width`](@ref) entries. Convenient at a call site that runs once. A particle loop should hold its own buffer and call -[`evaluate_all!`](@ref), which allocates nothing. +[`evaluate_all!`](@ref), which allocates nothing, or [`evaluate!`](@ref) for the value of a +spline at the point. """ function evaluate_all(b::AbstractBSplineBasis{T}, x::Number, d::Integer = 0) where {T} R = _evaltype(T, typeof(x)) diff --git a/src/spline.jl b/src/spline.jl index 4dc03cc..32cb224 100644 --- a/src/spline.jl +++ b/src/spline.jl @@ -119,6 +119,10 @@ The `d`-th derivative of `s` at `x`. For a tensor product `d` is a per-axis tupl evaluate(s::Spline, x) = evaluate(s.basis, s.coefficients, x) evaluate(s::Spline, x, d) = evaluate(s.basis, s.coefficients, x, d) +function evaluate!(buf::AbstractVector, s::Spline, x::Number, d::Integer = 0) + evaluate!(buf, s.basis, s.coefficients, x, d) +end + @doc raw""" SplineDerivative diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index b49565d..4d0d187 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -212,14 +212,36 @@ end function evaluate(B::TensorProductBasis{T, D}, û::AbstractArray{S, D}, x, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {T, D, S} + R = _evaltype(promote_type(T, S), eltype(x)) + bufs = ntuple(k -> Vector{R}(undef, local_width(B.bases[k])), D) + evaluate!(bufs, B, û, x, d) +end + +@doc raw""" + evaluate!(bufs::NTuple{D,AbstractVector}, B::TensorProductBasis, û::AbstractArray, x, + d = ntuple(_ -> 0, D)) + +The value of [`evaluate`](@ref)`(B, û, x, d)`, computed in the caller's per-axis buffers, +which allocates nothing. + +`bufs[k]` must have `local_width(bases(B)[k])` entries, or a `DimensionMismatch` is raised; +their contents on entry do not matter, and on return they are the blocks of +[`evaluate_all!`](@ref). Their element type is the type the sum is accumulated in, so it must +hold the promoted type of the basis, the coefficients and `x`. + +The buffers belong to one task. Two tasks that evaluate at once need a tuple each. +""" +function evaluate!(bufs::NTuple{D, <:AbstractVector{R}}, B::TensorProductBasis{T, D}, + û::AbstractArray{S, D}, x, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {R, T, D, S} size(û) == size(B) || throw(DimensionMismatch( "the coefficient array is $(size(û)) but the basis is $(size(B))")) - R = _evaltype(promote_type(T, S), eltype(x)) + widths = local_width(B) + map(length, bufs) == widths || throw(DimensionMismatch( + "evaluate! needs buffers of $(widths) entries on the axes, got $(map(length, bufs))")) # Only the local block contributes, so this is the fast path rather than a sum over the # whole array. `evaluate_all!` gives the p_d+1 nonzero factors on each axis; their outer # product against the corresponding block of û is the value. - bufs = ntuple(k -> Vector{R}(undef, local_width(B.bases[k])), D) j₀ = evaluate_all!(bufs, B, x, d) # `size(B, k)` reads `B.bases[k]`, and the axes need not have the same basis type, so @@ -399,52 +421,114 @@ mass_factors(op::KroneckerMass) = op.ops return before, after end -# Apply a one-dimensional operator along axis k of a D-dimensional array, in place. The array -# is viewed as (before, N_k, after) and each fibre along the middle index is solved on its own. +# Apply a one-dimensional operator along axis k of the D-dimensional array of size `dims` +# that `A` holds in its linear order, in place. The array is read as (before, N_k, after) and +# each fibre along the middle index is solved on its own. The linear index is computed rather +# than taken from a `reshape`, because a `reshape` of an `Array` allocates its header. # # The fibre is copied into a contiguous buffer rather than passed as a view. A view into the # middle index of a three-way reshape has stride `before`, and an FFTW plan or a LAPACK banded # solve addresses only contiguous memory -- so the one-dimensional solve would stage such a -# fibre through its own buffer in any case. The two buffers here are allocated once per axis -# and reused across that axis's fibres, so the cost is O(D) buffers per call rather than one -# pair per fibre. -# -# This adds no shared state of its own -- both buffers are local to the call -- but that is -# not thread safety: a `CirculantMass` holds its own scratch vector, which `mass_solve!` -# writes, so two threads applying the same operator would collide inside it. +# fibre through its own buffer in any case. The two buffers `src` and `dst` have N_k entries +# and are reused across the fibres of the axis. They can come from a caller's workspace, so +# the copies keep their bounds checks. # # `f!` receives (out, op, in) with two *distinct* buffers, so a method that cannot alias its # arguments -- `mul!` -- needs no defensive copy of its own. -function _apply_along!(f!, A::AbstractArray{T}, op, k::Integer) where {T} - dims = size(A) +function _apply_along!(f!, A::AbstractArray, dims::NTuple{D, Int}, op, k::Integer, + src::AbstractVector, dst::AbstractVector) where {D} before, after = _splitdims(dims, k) nk = dims[k] - A3 = reshape(A, before, nk, after) - src = Vector{T}(undef, nk) - dst = Vector{T}(undef, nk) + o = firstindex(A) - 1 for j in 1:after, i in 1:before - @inbounds for r in 1:nk - src[r] = A3[i, r, j] + base = o + i + before * nk * (j - 1) + for r in 1:nk + src[r] = A[base + before * (r - 1)] end f!(dst, op, src) - @inbounds for r in 1:nk - A3[i, r, j] = dst[r] + for r in 1:nk + A[base + before * (r - 1)] = dst[r] end end return A end -function mass_solve!(y::AbstractArray, op::KroneckerMass{T, D}, - x::AbstractArray) where {T, D} +# The form with buffers of its own, allocated once for the axis. +function _apply_along!(f!, A::AbstractArray{T}, op, k::Integer) where {T} + nk = size(A, k) + _apply_along!(f!, A, size(A), op, k, Vector{T}(undef, nk), Vector{T}(undef, nk)) +end + +# The solves along the axes, one axis per call. Recursion on the tuples rather than a loop +# over `1:D`: the axes need not have the same operator type, and indexing that tuple with a +# loop variable dispatches dynamically, which allocates. +_solve_axes!(y, dims, ::Tuple{}, ::Tuple{}, ::Tuple{}, k) = y +function _solve_axes!(y, dims, ops::Tuple, src::Tuple, dst::Tuple, k) + _apply_along!(mass_solve!, y, dims, first(ops), k, first(src), first(dst)) + _solve_axes!(y, dims, Base.tail(ops), Base.tail(src), Base.tail(dst), k + 1) +end + +struct KroneckerMassWorkspace{T, D} + src::NTuple{D, Vector{T}} + dst::NTuple{D, Vector{T}} +end + +Base.eltype(::KroneckerMassWorkspace{T}) where {T} = T + +@doc raw""" + mass_workspace(op::KroneckerMass, T) + mass_solve!(y, op::KroneckerMass, x, work) + +A workspace of element type `T` for [`mass_solve!`](@ref) on `op`, and the solve that uses it, +which allocates nothing. + +The workspace holds two buffers of `size(mass_factors(op)[k], 1)` entries for each axis `k`. +The solve takes `y` and `x` as arrays of the size of the tensor-product basis, or as vectors of +`size(op, 1)` entries, and `y === x` is allowed. A `y` of another length raises a +`DimensionMismatch`; a workspace whose element type is not `eltype(y)`, or whose buffer +lengths do not fit `op`, raises an `ArgumentError`. `mass_solve!(y, op, x)` builds a +workspace of `eltype(y)` on each call. + +A workspace belongs to one task. It does not make the operator safe for two tasks: a +[`CirculantMass`](@ref) factor holds a scratch vector of its own, which every solve writes. + +```jldoctest +julia> B = BSplineBasis(UniformMesh(8, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(6, 0 .. 1), 2); + +julia> op = mass_operator(TensorProductQuadrature(B)); x = randn(size(B)...); + +julia> work = mass_workspace(op, Float64); + +julia> mass_solve!(similar(x), op, x, work) == op \ x +true +``` +""" +function mass_workspace(op::KroneckerMass{T, D}, ::Type{S}) where {T, D, S} + KroneckerMassWorkspace{S, D}(map(n -> Vector{S}(undef, n), op.dims), + map(n -> Vector{S}(undef, n), op.dims)) +end + +function mass_solve!(y::AbstractArray, op::KroneckerMass{T, D}, x::AbstractArray, + work::KroneckerMassWorkspace) where {T, D} + eltype(work) === eltype(y) || throw(ArgumentError( + "the workspace has element type $(eltype(work)) but the array has element type " * + "$(eltype(y))")) + (map(length, work.src) == op.dims && map(length, work.dst) == op.dims) || + throw(ArgumentError( + "the workspace has buffers of lengths $(map(length, work.src)) but the operator " * + "needs $(op.dims)")) + length(y) == prod(op.dims) || throw(DimensionMismatch( + "the array has $(length(y)) entries but the operator acts on $(op.dims)")) y === x || copyto!(y, x) - Y = reshape(y, op.dims) - for k in 1:D - _apply_along!(mass_solve!, Y, op.ops[k], k) - end + _solve_axes!(y, op.dims, op.ops, work.src, work.dst, 1) return y end +function mass_solve!(y::AbstractArray, op::KroneckerMass, x::AbstractArray) + mass_solve!(y, op, x, mass_workspace(op, eltype(y))) +end + function LinearAlgebra.mul!(y::AbstractArray, op::KroneckerMass{T, D}, x::AbstractArray) where {T, D} y === x || copyto!(y, x) @@ -654,44 +738,155 @@ the ``O(N \prod_k n_k n_{q,k})`` of forming each entry by its own quadrature loo """ function contract(q::TensorProductQuadrature{T, D}, F::AbstractArray{S, D}, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {T, D, S} - size(F) == quadrature_grid_size(q) || throw(DimensionMismatch( - "the sample is $(size(F)) but the quadrature grid is $(quadrature_grid_size(q))")) - R = promote_type(T, S) + L = Array{R}(undef, size(q.basis)) + contract!(L, q, F, contract_workspace(q, R), d) +end + +struct ContractWorkspace{T, MW <: KroneckerMassWorkspace{T}} + a::Vector{T} + b::Vector{T} + mass::MW +end + +Base.eltype(::ContractWorkspace{T}) where {T} = T + +# The largest intermediate array of the contraction. Step k reads an array of +# N_1 ⋯ N_(k-1) Q_k ⋯ Q_D entries and writes one of N_1 ⋯ N_k Q_(k+1) ⋯ Q_D. +function _contract_length(q::TensorProductQuadrature{T, D}) where {T, D} + Qs = quadrature_grid_size(q) + Ns = size(q.basis) + n = m = prod(Qs) + for k in 1:D + m = m ÷ Qs[k] * Ns[k] + n = max(n, m) + end + return n +end + +""" + contract_workspace(q::TensorProductQuadrature, T) - # A copy, always. `convert(Array{R}, F)` is the *identity* when `F` is already an - # `Array{R}`, and the weighting below is in place -- so converting here would scale the - # caller's own sample by the quadrature weights and leave it that way. `F` is very often - # exactly such an array, since `quadrature_sample` returns one. - B = Array{R}(undef, size(F)) - copyto!(B, F) +A workspace of element type `T` for [`contract!`](@ref) and the four-argument +[`l2_projection!`](@ref) on `q`. A workspace belongs to one task. +""" +function contract_workspace(q::TensorProductQuadrature, ::Type{S}) where {S} + n = _contract_length(q) + ContractWorkspace(Vector{S}(undef, n), Vector{S}(undef, n), + mass_workspace(q.mass, S)) +end + +@doc raw""" + contract!(L, q::TensorProductQuadrature, F, work, d = ntuple(_ -> 0, D)) + l2_projection!(û, q::TensorProductQuadrature, F::AbstractArray, work) + +The forms of [`contract`](@ref) and [`l2_projection!`](@ref) that take a workspace from +[`contract_workspace`](@ref), and allocate nothing. + +`contract!` writes [`contract`](@ref)`(q, F, d)` into `L`, of size `size(basis(q))`; +`l2_projection!` writes `l2_projection(q, F)` into `û`. `F` is the sample on the quadrature +grid and is not changed. The workspace holds two buffers for the intermediate arrays of the +contraction and a [`mass_workspace`](@ref) of `mass_operator(q)` for the solve. Wrong sizes +of `L`, `û` or `F` raise the `DimensionMismatch` of the allocating forms; a workspace whose +element type is not `eltype(L)`, or `eltype(û)`, or whose sizes do not fit `q`, raises an +`ArgumentError`. `contract` builds a workspace on each call, and so does the three-argument +`l2_projection!`. + +A workspace belongs to one task. It does not make the quadrature safe for two tasks: a +[`CirculantMass`](@ref) factor of its mass operator holds a scratch vector of its own. + +```jldoctest +julia> B = BSplineBasis(UniformMesh(8, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(6, 0 .. 1), 2); + +julia> q = TensorProductQuadrature(B); F = quadrature_sample(q, x -> x[1] * x[2]); + +julia> work = contract_workspace(q, Float64); L = zeros(size(B)...); + +julia> contract!(L, q, F, work) == contract(q, F) +true +``` +""" +function contract!(L::AbstractArray, q::TensorProductQuadrature{T, D}, + F::AbstractArray{S, D}, work::ContractWorkspace, + d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {T, D, S} + size(F) == quadrature_grid_size(q) || throw(DimensionMismatch( + "the sample is $(size(F)) but the quadrature grid is $(quadrature_grid_size(q))")) + size(L) == size(q.basis) || throw(DimensionMismatch( + "the load array is $(size(L)) but the basis is $(size(q.basis))")) + eltype(work) === eltype(L) || throw(ArgumentError( + "the workspace has element type $(eltype(work)) but the load array has element " * + "type $(eltype(L))")) + n = _contract_length(q) + (length(work.a) == n && length(work.b) == n) || throw(ArgumentError( + "the workspace holds buffers of $(length(work.a)) entries but this quadrature " * + "needs $(n)")) + + # A copy, always: the weighting below is in place, and `F` is the caller's sample. + copyto!(work.a, firstindex(work.a), F, firstindex(F), length(F)) # The weights of every axis, applied once, before any contraction. Doing it here rather # than inside the loop keeps the weighting a single elementwise pass and makes the # per-axis step a plain matrix product. - for k in 1:D - w = quadrature_weights(q.quadratures[k]) - _scale_along!(B, w, k) - end + _scale_axes!(work.a, size(F), q.quadratures, 1) - for k in 1:D - Φ = basis_values(q.quadratures[k], d[k]) - Q = size(Φ, 2) - M = reshape(B, Q, :) - C = Φ * M # (N_k, rest) - B = permutedims(C) # (rest, N_k) -- cycles this axis to the back - end + C = _contract_axes!(work.a, work.b, length(F), q.quadratures, d) + copyto!(L, firstindex(L), C, firstindex(C), length(L)) + return L +end - return reshape(B, size(q.basis)) +# The weighting and the contraction, one axis per call. Recursion on the tuple rather than a +# loop over `1:D`, for the reason of `_solve_axes!`: the axes need not have the same +# quadrature type. +_scale_axes!(A, dims, ::Tuple{}, k) = A +function _scale_axes!(A, dims, qs::Tuple, k) + _scale_along!(A, dims, quadrature_weights(first(qs)), k) + _scale_axes!(A, dims, Base.tail(qs), k + 1) end -function _scale_along!(A::AbstractArray, w::AbstractVector, k::Integer) - dims = size(A) +# `src` holds an array of `n` entries whose first axis is the quadrature axis of `first(qs)`. +# Each step writes `Φ * M` into `dst` transposed, as a (rest, N_k) array, which cycles the +# contracted axis to the back; after the last step the axes are back in order, and the +# buffer that holds the result is returned. +_contract_axes!(src, dst, n, ::Tuple{}, ::Tuple{}) = src +function _contract_axes!(src, dst, n, qs::Tuple, d::Tuple) + Φ = basis_values(first(qs), first(d)) + N, Q = size(Φ) + rest = n ÷ Q + _contract_axis!(dst, Φ, src, rest) + _contract_axes!(dst, src, rest * N, Base.tail(qs), Base.tail(d)) +end + +# dst[r, i] = Σ_q Φ[i, q] src[q, r], with `src` read as (Q, rest) and `dst` written as +# (rest, N). The sum runs over the stored entries of each column of Φ, with `muladd`, in the +# order of SparseArrays' sparse-times-dense product, so the result is the one that +# `permutedims(Φ * M)` gives. +function _contract_axis!(dst::AbstractVector, Φ::SparseMatrixCSC, src::AbstractVector, + rest::Integer) + N, Q = size(Φ) + rows = rowvals(Φ) + vals = nonzeros(Φ) + @inbounds for r in 1:rest + for i in 1:N + dst[r + rest * (i - 1)] = zero(eltype(dst)) + end + for c in 1:Q + s = src[c + Q * (r - 1)] + for t in nzrange(Φ, c) + i = r + rest * (rows[t] - 1) + dst[i] = muladd(vals[t], s, dst[i]) + end + end + end + return dst +end + +# Scale axis k of the array of size `dims` that `A` holds in its linear order by `w`. +function _scale_along!(A::AbstractVector, dims::NTuple{D, Int}, w::AbstractVector, + k::Integer) where {D} before, after = _splitdims(dims, k) nk = dims[k] - A3 = reshape(A, before, nk, after) @inbounds for j in 1:after, r in 1:nk, i in 1:before - A3[i, r, j] *= w[r] + A[i + before * ((r - 1) + nk * (j - 1))] *= w[r] end return A end @@ -699,12 +894,15 @@ end @doc raw""" l2_projection(q::TensorProductQuadrature, f) l2_projection!(û, q::TensorProductQuadrature, f) + l2_projection!(û, q::TensorProductQuadrature, F::AbstractArray, work) The coefficient array of the ``L^2`` projection of `f` onto the tensor-product spline space, ``\hat{u} = \mathbb{M}^{-1} \int_\Omega f \, \Phi_I \, dx``. `f` may be a function of a `D`-tuple of coordinates, or an array already sampled on the -quadrature grid — see [`quadrature_sample`](@ref). The result has size `size(basis(q))`. +quadrature grid — see [`quadrature_sample`](@ref). The result has size `size(basis(q))`. The +form with a workspace from [`contract_workspace`](@ref) takes the sample and allocates +nothing. ```jldoctest julia> B = BSplineBasis(UniformMesh(16, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(16, 0 .. 1), 3); @@ -732,10 +930,15 @@ function l2_projection(q::TensorProductQuadrature{T, D}, f) where {T, D} end function l2_projection!(û::AbstractArray, q::TensorProductQuadrature{T, D}, f) where {T, D} + F = f isa AbstractArray ? f : quadrature_sample(q, f) + l2_projection!(û, q, F, contract_workspace(q, eltype(û))) +end + +function l2_projection!(û::AbstractArray, q::TensorProductQuadrature{T, D}, + F::AbstractArray, work::ContractWorkspace) where {T, D} size(û) == size(q.basis) || throw(DimensionMismatch( "the coefficient array is $(size(û)) but the basis is $(size(q.basis))")) - F = f isa AbstractArray ? f : quadrature_sample(q, f) - copyto!(û, contract(q, F)) - mass_solve!(û, mass_operator(q), û) + contract!(û, q, F, work) + mass_solve!(û, mass_operator(q), û, work.mass) return û end diff --git a/test/Project.toml b/test/Project.toml index 643716d..a291927 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,7 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -13,6 +14,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" Documenter = "1" +ForwardDiff = "1" SafeTestsets = "0.1" StaticArrays = "1" Test = "1" diff --git a/test/integration/workspaces.jl b/test/integration/workspaces.jl new file mode 100644 index 0000000..8abff93 --- /dev/null +++ b/test/integration/workspaces.jl @@ -0,0 +1,347 @@ +using ForwardDiff: Dual +using LinearAlgebra +using Random +using SimpleSplines +using Test + +Random.seed!(1234) + +# Fixed-arity barriers: the arguments arrive with concrete types, and the call runs once +# before it is measured. +allocs(f::F, a::A, b::B, c::C) where {F, A, B, C} = (f(a, b, c); @allocated f(a, b, c)) +function allocs(f::F, a::A, b::B, c::C, d::D) where {F, A, B, C, D} + f(a, b, c, d) + @allocated f(a, b, c, d) +end +function allocs(f::F, a::A, b::B, c::C, d::D, e::E) where {F, A, B, C, D, E} + f(a, b, c, d, e) + @allocated f(a, b, c, d, e) +end + +# Coverage instrumentation allocates on every line, so the readings mean nothing under it. +const MEASURE = Base.JLOptions().code_coverage == 0 + +# The outcome of a call: its value, or the type of the exception it raised. +outcome(f) = + try + f() + catch e + typeof(e) + end + +# The message of the exception that `f` raises. +message(f) = + try + f() + "" + catch e + sprint(showerror, e) + end + +const Dual12 = Dual{Nothing, Float64, 12} +dual(x) = Dual{Nothing}(Float64(x), ntuple(i -> Float64(i) / 7, 12)...) + +function bases_1d(::Type{T}) where {T} + m = UniformMesh(8, T(0) .. T(1)) + (BSplineBasis(m, 1), BSplineBasis(m, 3), BSplineBasis(m, 3, Periodic()), + BSplineBasis(m, 3, Dirichlet()), BSplineBasis(m, 3, Natural())) +end + +function basis_2d(::Type{T}, n = 8) where {T} + BSplineBasis(UniformMesh(n, T(0) .. T(1)), 3, Periodic()) ⊗ + BSplineBasis(UniformMesh(n, T(-1) .. T(2)), 3) +end + +function basis_3d(::Type{T}, n = 4) where {T} + BSplineBasis(UniformMesh(n, T(0) .. T(1)), 3) ⊗ + BSplineBasis(UniformMesh(n, T(0) .. T(2)), 3, Periodic()) ⊗ + BSplineBasis(UniformMesh(n, T(-1) .. T(1)), 3, Dirichlet()) +end + +# The same basis on more cells: a workspace built for it has longer buffers, so a missing +# size check shows as a wrong result or an error rather than as a write out of bounds. +larger(B) = ndims(B) == 2 ? basis_2d(eltype(B), 9) : basis_3d(eltype(B), 5) + +@testset "$(rpad("evaluate! on a one-dimensional basis",76))" begin + for T in (Float32, Float64), b in bases_1d(T) + + û = randn(T, nbasis(b)) + buf = zeros(T, local_width(b)) + s = Spline(b, û) + # interior points, a breakpoint, both ends, and points outside the domain + for x in T[0.37, 0.5, 0, 1, prevfloat(T(1)), -0.25, 1.75], d in 0:(degree(b) + 1) + + @test (@inferred evaluate!(buf, b, û, x, d)) == evaluate(b, û, x, d) + @test (@inferred evaluate!(buf, s, x, d)) == evaluate(s, x, d) + end + @test evaluate!(buf, b, û, T(0.37)) == evaluate(b, û, T(0.37)) + @test evaluate!(buf, s, T(0.37)) == s(T(0.37)) + + # A Dual argument with a Dual buffer. `isequal` rather than `==`: a periodic basis at + # the ends of its domain gives NaN partials on both sides (K6 of KNOWN_ISSUES.md), and + # `==` is false for them. + bufd = zeros(Dual12, local_width(b)) + for x in (0.37, 0.5, 1.0, 1.75), d in 0:1 + + @test isequal((@inferred evaluate!(bufd, b, û, dual(x), d)), + evaluate(b, û, dual(x), d)) + @test isequal((@inferred evaluate!(bufd, s, dual(x), d)), evaluate(s, dual(x), d)) + end + end +end + +@testset "$(rpad("evaluate! allocates nothing",76))" begin + if MEASURE + for T in (Float32, Float64), b in bases_1d(T), d in 0:1 + û = randn(T, nbasis(b)) + s = Spline(b, û) + x = T(0.37) + buf = zeros(T, local_width(b)) + @test allocs(evaluate!, buf, b, û, x, d) == 0 + @test allocs(evaluate!, buf, s, x, d) == 0 + @test allocs(evaluate, b, û, x, d) > 0 + @test allocs(evaluate, s, x, d) > 0 + + bufd = zeros(Dual12, local_width(b)) + @test allocs(evaluate!, bufd, b, û, dual(0.37), d) == 0 + @test allocs(evaluate!, bufd, s, dual(0.37), d) == 0 + @test allocs(evaluate, b, û, dual(0.37), d) > 0 + end + end +end + +@testset "$(rpad("evaluate! refuses a buffer of the wrong length",76))" begin + for b in bases_1d(Float64) + û = randn(nbasis(b)) + w = local_width(b) + for n in (w - 1, w + 1) + @test_throws DimensionMismatch evaluate!(zeros(n), b, û, 0.37) + msg = message(() -> evaluate!(zeros(n), b, û, 0.37)) + @test startswith(msg, + "DimensionMismatch: evaluate! needs a buffer of $(w) entries for this " * + "basis, got $(n)") + @test_throws DimensionMismatch evaluate!(zeros(n), Spline(b, û), 0.37) + end + # the coefficient check of `evaluate`, unchanged + @test_throws DimensionMismatch evaluate!(zeros(w), b, zeros(nbasis(b) + 1), 0.37) + @test startswith(message(() -> evaluate!(zeros(w), b, zeros(nbasis(b) + 1), 0.37)), + "DimensionMismatch: the coefficient vector has $(nbasis(b) + 1) entries") + end +end + +@testset "$(rpad("evaluate! on hostile arguments matches evaluate",76))" begin + for T in (Float32, Float64), b in bases_1d(T) + + û = randn(T, nbasis(b)) + buf = zeros(T, local_width(b)) + for x in T[0, -0.0, nextfloat(T(0)), floatmin(T), eps(T), 1, floatmax(T) / 2, + floatmax(T), Inf, -Inf, NaN, 1e-8, 1e-4, 1e4, 1e8], + d in (0, 1, 4, -1) + + @test isequal(outcome(() -> evaluate!(buf, b, û, x, d)), + outcome(() -> evaluate(b, û, x, d))) + end + # `d < 0` is the `ArgumentError` of `evaluate_all!` + @test_throws ArgumentError evaluate!(buf, b, û, T(0.5), -1) + end + # BigFloat: the values agree + for b in bases_1d(BigFloat) + û = randn(nbasis(b)) + buf = zeros(BigFloat, local_width(b)) + for x in big.([0.37, 1.0, 1.5]), d in 0:2 + + @test evaluate!(buf, b, û, x, d) == evaluate(b, û, x, d) + end + end + # degree 0 on one cell + for b in (BSplineBasis(UniformMesh(1, 0 .. 1), 0), + BSplineBasis(UniformMesh(1, 0 .. 1), 0, Periodic())) + û = [2.5] + for x in (0.0, 0.5, 1.0, 1.5, -0.5), d in 0:1 + + @test evaluate!(zeros(1), b, û, x, d) == evaluate(b, û, x, d) + end + end +end + +@testset "$(rpad("evaluate! on a tensor-product basis",76))" begin + for T in (Float32, Float64), B in (basis_2d(T), basis_3d(T)) + + D = ndims(B) + û = randn(T, size(B)...) + bufs = map(w -> zeros(T, w), local_width(B)) + bufsd = map(w -> zeros(Dual12, w), local_width(B)) + points = D == 2 ? (T.((0.37, 0.4)), T.((0.0, 2.0)), T.((1.5, 3.0))) : + (T.((0.37, 1.2, -0.3)), T.((1.0, 0.0, 1.0)), T.((-0.5, 2.5, 0.2))) + orders = (ntuple(_ -> 0, D), ntuple(k -> k == 1 ? 1 : 0, D), + ntuple(k -> k == D ? 1 : 0, D)) + for x in points, d in orders + + @test (@inferred evaluate!(bufs, B, û, x, d)) == evaluate(B, û, x, d) + xd = map(dual, x) + @test isequal((@inferred evaluate!(bufsd, B, û, xd, d)), evaluate(B, û, xd, d)) + end + @test evaluate!(bufs, B, û, first(points)) == evaluate(B, û, first(points)) + + if MEASURE + x = first(points) + d = ntuple(k -> k == 1 ? 1 : 0, D) + @test allocs(evaluate!, bufs, B, û, x, d) == 0 + @test allocs(evaluate, B, û, x, d) > 0 + @test allocs(evaluate!, bufsd, B, û, map(dual, x), d) == 0 + @test allocs(evaluate, B, û, map(dual, x), d) > 0 + end + + # a wrong length on any one axis + for k in 1:D + w = local_width(B) + bad = ntuple(i -> zeros(T, i == k ? w[i] + 1 : w[i]), D) + @test_throws DimensionMismatch evaluate!(bad, B, û, first(points)) + @test startswith(message(() -> evaluate!(bad, B, û, first(points))), + "DimensionMismatch: evaluate! needs buffers of $(w) entries on the axes, " * + "got $(map(length, bad))") + end + # a coefficient array of the wrong size, as `evaluate` refuses it + ûbad = zeros(T, (size(B) .+ 1)...) + @test_throws DimensionMismatch evaluate!(bufs, B, ûbad, first(points)) + @test startswith(message(() -> evaluate!(bufs, B, ûbad, first(points))), + "DimensionMismatch: the coefficient array is $(size(ûbad)) but the basis is " * + "$(size(B))") + end +end + +@testset "$(rpad("mass_solve! with a workspace",76))" begin + for T in (Float32, Float64), B in (basis_2d(T), basis_3d(T)) + + op = mass_operator(TensorProductQuadrature(B)) + x = randn(T, size(B)...) + work = mass_workspace(op, T) + y = similar(x) + @test (@inferred mass_solve!(y, op, x, work)) === y + @test y == mass_solve!(similar(x), op, x) + # an independent reference: the assembled Kronecker product + @test vec(y) ≈ Matrix(op) \ vec(x) rtol = 100 * eps(T) * cond(Matrix(op)) + # y and x may alias, and a vector of length `length(B)` is accepted + z = copy(x) + @test mass_solve!(z, op, z, work) == y + v = vec(copy(x)) + @test mass_solve!(v, op, v, work) == vec(y) + + if MEASURE + @test allocs(mass_solve!, y, op, x, work) == 0 + @test allocs(mass_solve!, y, op, x) > 0 + end + + # a workspace of the wrong element type names both types + S = T === Float64 ? Float32 : Float64 + bad = mass_workspace(op, S) + @test_throws ArgumentError mass_solve!(y, op, x, bad) + @test startswith(message(() -> mass_solve!(y, op, x, bad)), + "ArgumentError: the workspace has element type $(S) but the array has element " * + "type $(T)") + # a workspace of a larger operator names both sizes + other = mass_operator(TensorProductQuadrature(larger(B))) + badsize = mass_workspace(other, T) + @test_throws ArgumentError mass_solve!(y, op, x, badsize) + @test startswith(message(() -> mass_solve!(y, op, x, badsize)), + "ArgumentError: the workspace has buffers of lengths " * + "$(map(n -> size(n, 1), mass_factors(other))) but the operator needs " * + "$(map(n -> size(n, 1), mass_factors(op)))") + # a workspace of another dimension + # (of 19 functions per axis, longer than every axis of `B`) + b19 = BSplineBasis(UniformMesh(16, T(0) .. T(1)), 3) + flat = mass_workspace(mass_operator(TensorProductQuadrature(b19 ⊗ b19 ⊗ b19 ⊗ b19)), T) + @test_throws ArgumentError mass_solve!(y, op, x, flat) + # a long `dst` buffer + longdst = SimpleSplines.KroneckerMassWorkspace(work.src, + map(v -> [v; zero(T)], work.dst)) + @test_throws ArgumentError mass_solve!(y, op, x, longdst) + @test startswith(message(() -> mass_solve!(y, op, x, longdst)), + "ArgumentError: the workspace has buffers of lengths") + # an array of the wrong length + long = zeros(T, length(x) + 1) + @test_throws DimensionMismatch mass_solve!(long, op, vec(x), work) + @test startswith(message(() -> mass_solve!(long, op, vec(x), work)), + "DimensionMismatch: the array has $(length(x) + 1) entries but the operator " * + "acts on $(size(x))") + end +end + +# The load array by its definition, from the dense per-axis tables. +function reference_contract(q, F, d) + A = [Matrix(basis_values(quadratures(q)[k], d[k])) * + Diagonal(quadrature_weights(quadratures(q)[k])) for k in 1:ndims(q)] + reshape(foldl(kron, reverse(A)) * vec(F), size(q)...) +end + +@testset "$(rpad("contract! and l2_projection! with a workspace",76))" begin + for T in (Float32, Float64), B in (basis_2d(T), basis_3d(T)) + + D = ndims(B) + q = TensorProductQuadrature(B) + F = randn(T, quadrature_grid_size(q)...) + F₀ = copy(F) + work = contract_workspace(q, T) + L = zeros(T, size(B)...) + for d in (ntuple(_ -> 0, D), ntuple(k -> k == D ? 1 : 0, D)) + @test (@inferred contract!(L, q, F, work, d)) === L + @test L == contract(q, F, d) + @test L ≈ reference_contract(q, F, d) rtol = 100 * eps(T) + end + @test contract!(L, q, F, work) == contract(q, F) + @test F == F₀ # the sample is not consumed + + û = zeros(T, size(B)...) + @test (@inferred l2_projection!(û, q, F, work)) === û + @test û == l2_projection!(zeros(T, size(B)...), q, F) + @test û == l2_projection(q, F) + @test F == F₀ + + if MEASURE + d = ntuple(k -> k == D ? 1 : 0, D) + @test allocs(contract!, L, q, F, work, d) == 0 + @test allocs(contract, q, F, d) > 0 + @test allocs(l2_projection!, û, q, F, work) == 0 + @test allocs(l2_projection!, û, q, F) > 0 + end + + # wrong sizes, as the allocating forms raise them + Fbad = zeros(T, (quadrature_grid_size(q) .+ 1)...) + @test_throws DimensionMismatch contract!(L, q, Fbad, work) + @test message(() -> contract!(L, q, Fbad, work)) == + message(() -> contract(q, Fbad)) + Lbad = zeros(T, (size(B) .+ 1)...) + @test_throws DimensionMismatch contract!(Lbad, q, F, work) + @test startswith(message(() -> contract!(Lbad, q, F, work)), + "DimensionMismatch: the load array is $(size(Lbad)) but the basis is $(size(B))") + @test_throws DimensionMismatch l2_projection!(Lbad, q, F, work) + @test startswith(message(() -> l2_projection!(Lbad, q, F, work)), + "DimensionMismatch: the coefficient array is $(size(Lbad)) but the basis is " * + "$(size(B))") + + # a workspace of the wrong element type names both types + S = T === Float64 ? Float32 : Float64 + bad = contract_workspace(q, S) + @test_throws ArgumentError contract!(L, q, F, bad) + @test startswith(message(() -> contract!(L, q, F, bad)), + "ArgumentError: the workspace has element type $(S) but the load array has " * + "element type $(T)") + @test_throws ArgumentError l2_projection!(û, q, F, bad) + # a workspace of a larger quadrature names both lengths + badsize = contract_workspace(TensorProductQuadrature(larger(B)), T) + @test_throws ArgumentError contract!(L, q, F, badsize) + @test startswith(message(() -> contract!(L, q, F, badsize)), + "ArgumentError: the workspace holds buffers of $(length(badsize.a)) entries but " * + "this quadrature needs $(length(work.a))") + @test length(badsize.a) != length(work.a) + @test_throws ArgumentError l2_projection!(û, q, F, badsize) + # a long second buffer + longb = SimpleSplines.ContractWorkspace(work.a, [work.b; zero(T)], work.mass) + @test_throws ArgumentError contract!(L, q, F, longb) + # a mass workspace of a larger operator + wrongmass = SimpleSplines.ContractWorkspace(work.a, work.b, badsize.mass) + @test_throws ArgumentError l2_projection!(û, q, F, wrongmass) + @test startswith(message(() -> l2_projection!(û, q, F, wrongmass)), + "ArgumentError: the workspace has buffers of lengths") + end +end diff --git a/test/quality/jet.jl b/test/quality/jet.jl index 81f87da..510e33d 100644 --- a/test/quality/jet.jl +++ b/test/quality/jet.jl @@ -70,6 +70,31 @@ if isdefined(JET, :JET_AVAILABLE) ? JET.JET_AVAILABLE : JET.JET_LOADABLE (S2, typeof(opb), Vector{Float64}); target_modules = (SimpleSplines,)))) @test isempty(JET.get_reports(JET.report_opt(mass_solve!, (S2, typeof(opc), Vector{Float64}); target_modules = (SimpleSplines,)))) + + # test/integration/workspaces.jl: the caller-held buffers and workspaces, on a periodic × + # clamped product, whose axes differ in type + m8 = UniformMesh(8, 0 .. 1) + for b1 in (BSplineBasis(m8, 3), BSplineBasis(m8, 3, Periodic()), + BSplineBasis(m8, 3, Natural())) + @test isempty(JET.get_reports(JET.report_opt(evaluate!, + (Vector{Float64}, typeof(b1), Vector{Float64}, Float64, Int); + target_modules = (SimpleSplines,)))) + end + Bw = BSplineBasis(m8, 3, Periodic()) ⊗ BSplineBasis(UniformMesh(8, -1 .. 2), 3) + qw = TensorProductQuadrature(Bw) + @test isempty(JET.get_reports(JET.report_opt(evaluate!, + (NTuple{2, Vector{Float64}}, typeof(Bw), Matrix{Float64}, NTuple{2, Float64}, + NTuple{2, Int}); target_modules = (SimpleSplines,)))) + workw = contract_workspace(qw, Float64) + @test isempty(JET.get_reports(JET.report_opt(mass_solve!, + (Matrix{Float64}, typeof(mass_operator(qw)), Matrix{Float64}, typeof(workw.mass)); + target_modules = (SimpleSplines,)))) + @test isempty(JET.get_reports(JET.report_opt(contract!, + (Matrix{Float64}, typeof(qw), Matrix{Float64}, typeof(workw), NTuple{2, Int}); + target_modules = (SimpleSplines,)))) + @test isempty(JET.get_reports(JET.report_opt(l2_projection!, + (Matrix{Float64}, typeof(qw), Matrix{Float64}, typeof(workw)); + target_modules = (SimpleSplines,)))) else @test_skip "JET does not work on this Julia version" # aviatesk/JET.jl#681 end diff --git a/test/runtests.jl b/test/runtests.jl index 4396cfc..1cd9419 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,6 +13,7 @@ if "core" in GROUPS @safetestset "Spline quadrature" include("quadrature.jl") @safetestset "Tensor products" include("tensorproduct.jl") @safetestset "Polar splines" include("polar.jl") + @safetestset "Caller-held workspaces" include("integration/workspaces.jl") end if "doctests" in GROUPS @safetestset "Doctests" include("quality/doctests.jl") From eac719b7ddf69515aac3dc23d12aaef56aaa0d87 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 03:25:17 +0100 Subject: [PATCH 02/12] Name both buffer lengths in the workspace size errors The size refusals of mass_solve! and contract! name the lengths of both buffers, and the tests assert each at its own place, for a long first, a long second and two long buffers. The three-argument l2_projection! checks the size of its array before it samples f. The allocating contract builds its two buffers only, without a mass workspace. The weighting in contract! loops over the weights, which have one type. CHANGELOG: the bytes of contract, and the bit-for-bit claim limited to equal element types. KNOWN_ISSUES: K7, K8, K9. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 21 ++++++++------ KNOWN_ISSUES.md | 33 ++++++++++++++++++++++ src/tensorproduct.jl | 49 +++++++++++++++++++-------------- test/integration/workspaces.jl | 50 +++++++++++++++++++++++++++++----- 4 files changed, 116 insertions(+), 37 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9f4f8ea..6f79243 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,17 +28,20 @@ cells. Each pair is `Float64` / `Float32`. | `evaluate(b, û, x, d)` | 96 / 80 | 96 / 80 | 0 / 0 | | `evaluate(B, û, x, d)` | 192 / 160 | 192 / 160 | 0 / 0 | | `mass_solve!(y, op, x)` | 1232 / 1104 | 544 / 416 | 0 / 0 | -| `contract(q, F, d)` | 24064 / 13536 | 30128 / 15328 | 0 / 0 | +| `contract(q, F, d)` | 24064 / 13536 | 29584 / 14912 | 0 / 0 | | `l2_projection!(û, q, F)` | 25296 / 14640 | 29344 / 14880 | 0 / 0 | -The solve along the axes of a `KroneckerMass` and the weighting and contraction along the axes -of `contract` are unrolled over the tuple of axes, so an operator or a quadrature whose axes have -different types no longer dispatches dynamically: that is the drop of the three-argument -`mass_solve!` on the periodic × clamped basis. The allocating `contract` and `l2_projection!` -allocate more than before, because they build the whole workspace: two buffers of the largest -intermediate array, and the buffers of the mass solve. `mass_solve!` on a `KroneckerMass` checks -the length of `y` before it copies `x`, and raises a `DimensionMismatch` for a `y` of the wrong -length. The results of the allocating forms are bit-for-bit those of `cb82e6f`. +The solve along the axes of a `KroneckerMass` and the contraction along the axes of `contract` +are unrolled over the tuple of axes, so an operator or a quadrature whose axes have different +types no longer dispatches dynamically: that is the drop of the three-argument `mass_solve!` on +the periodic × clamped basis. The allocating `contract` and `l2_projection!` +allocate more than before, because each builds two buffers of the largest intermediate array, +and `l2_projection!` also the buffers of the mass solve. `mass_solve!` on a `KroneckerMass` +checks the length of `y` before it copies `x`, and raises a `DimensionMismatch` for a `y` of the +wrong length. Where the element types of the arguments are equal, the results of the allocating +forms are bit-for-bit those of `cb82e6f`. The three-argument `l2_projection!` now contracts in +`eltype(û)`, not in the promoted type of the quadrature and the sample, so where these types +differ its result can differ from that of `cb82e6f` by the round-off of the narrower type. `ForwardDiff` joins the test dependencies: the tests evaluate with `ForwardDiff.Dual` arguments and buffers, which also allocate nothing. diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 3fda417..df27dc7 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -85,3 +85,36 @@ are right; only the derivatives are lost. - **kind:** defect - **found:** 2026-11-02 + +### K7 · `l2_projection!` into a `Float32` array fails on a `Float64` periodic quadrature + +- **location:** `src/tensorproduct.jl`, `l2_projection!`, and the FFT solve of `CirculantMass` +- **evidence:** with `B = BSplineBasis(UniformMesh(8, 0.0 .. 1.0), 3, Periodic()) ⊗ + BSplineBasis(UniformMesh(8, -1.0 .. 2.0), 3)` and `q = TensorProductQuadrature(B)`, + `l2_projection!(zeros(Float32, size(B)...), q, ones(quadrature_grid_size(q)...))` raises + `MethodError: no method matching mul!(::Vector{ComplexF64}, ::FFTW.rFFTWPlan{Float64, -1, + false, 1, Tuple{Int64}}, ::Vector{Float32}, ::Bool, …`, on `cb82e6f` and on the branch that + adds caller-held workspaces alike: the `Float64` plan of the periodic factor does not take a + `Float32` vector. +- **kind:** defect +- **found:** 2026-11-02; older than `cb82e6f` + +### K8 · `mass_solve!` on a `KroneckerMass` does not check the length of `x` + +- **location:** `src/tensorproduct.jl`, `mass_solve!(y, op::KroneckerMass, x)` +- **evidence:** the solve copies `x` into `y` with `copyto!`, which accepts a shorter `x`, and + then solves on the whole of `y`. On the basis of K7, with `y = zeros(size(B)...)` (88 entries) + and `x = ones(83)`, `mass_solve!(y, mass_operator(q), x)` returns with no error, on `cb82e6f` + and on the branch that adds caller-held workspaces alike; the last five entries of `x` are + taken from the zeros of `y`. +- **kind:** defect +- **found:** 2026-11-02; older than `cb82e6f` + +### K9 · `evaluate!` on a `TensorProductBasis` takes buffers of one type only + +- **location:** `src/tensorproduct.jl`, `evaluate!(bufs, B::TensorProductBasis, û, x, d)` +- **evidence:** the method takes `bufs` as an `NTuple{D, <:AbstractVector{R}}`, so a tuple + of a `Vector` and a `SubArray`, or of vectors of different element types, raises a + `MethodError`. `evaluate_all!(bufs, B, x, d)` on `cb82e6f` has the same signature. +- **kind:** found late +- **found:** 2026-11-02 diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 4d0d187..42734ca 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -516,8 +516,8 @@ function mass_solve!(y::AbstractArray, op::KroneckerMass{T, D}, x::AbstractArray "$(eltype(y))")) (map(length, work.src) == op.dims && map(length, work.dst) == op.dims) || throw(ArgumentError( - "the workspace has buffers of lengths $(map(length, work.src)) but the operator " * - "needs $(op.dims)")) + "the workspace has buffers of lengths $(map(length, work.src)) and " * + "$(map(length, work.dst)) but the operator needs $(op.dims)")) length(y) == prod(op.dims) || throw(DimensionMismatch( "the array has $(length(y)) entries but the operator acts on $(op.dims)")) y === x || copyto!(y, x) @@ -740,10 +740,13 @@ function contract(q::TensorProductQuadrature{T, D}, F::AbstractArray{S, D}, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {T, D, S} R = promote_type(T, S) L = Array{R}(undef, size(q.basis)) - contract!(L, q, F, contract_workspace(q, R), d) + n = _contract_length(q) + # No mass workspace: `contract!` does not read it. + contract!( + L, q, F, ContractWorkspace(Vector{R}(undef, n), Vector{R}(undef, n), nothing), d) end -struct ContractWorkspace{T, MW <: KroneckerMassWorkspace{T}} +struct ContractWorkspace{T, MW <: Union{KroneckerMassWorkspace{T}, Nothing}} a::Vector{T} b::Vector{T} mass::MW @@ -789,8 +792,8 @@ grid and is not changed. The workspace holds two buffers for the intermediate ar contraction and a [`mass_workspace`](@ref) of `mass_operator(q)` for the solve. Wrong sizes of `L`, `û` or `F` raise the `DimensionMismatch` of the allocating forms; a workspace whose element type is not `eltype(L)`, or `eltype(û)`, or whose sizes do not fit `q`, raises an -`ArgumentError`. `contract` builds a workspace on each call, and so does the three-argument -`l2_projection!`. +`ArgumentError`. `contract` builds the two buffers on each call, and the three-argument +`l2_projection!` a whole workspace. A workspace belongs to one task. It does not make the quadrature safe for two tasks: a [`CirculantMass`](@ref) factor of its mass operator holds a scratch vector of its own. @@ -818,31 +821,28 @@ function contract!(L::AbstractArray, q::TensorProductQuadrature{T, D}, "type $(eltype(L))")) n = _contract_length(q) (length(work.a) == n && length(work.b) == n) || throw(ArgumentError( - "the workspace holds buffers of $(length(work.a)) entries but this quadrature " * - "needs $(n)")) + "the workspace holds buffers of $(length(work.a)) and $(length(work.b)) entries but " * + "this quadrature needs $(n)")) # A copy, always: the weighting below is in place, and `F` is the caller's sample. copyto!(work.a, firstindex(work.a), F, firstindex(F), length(F)) # The weights of every axis, applied once, before any contraction. Doing it here rather # than inside the loop keeps the weighting a single elementwise pass and makes the - # per-axis step a plain matrix product. - _scale_axes!(work.a, size(F), q.quadratures, 1) + # per-axis step a plain matrix product. The weights of all axes are `Vector{T}`, so the + # loop indexes a tuple of one type. + w = quadrature_weights(q) + for k in 1:D + _scale_along!(work.a, size(F), w[k], k) + end C = _contract_axes!(work.a, work.b, length(F), q.quadratures, d) copyto!(L, firstindex(L), C, firstindex(C), length(L)) return L end -# The weighting and the contraction, one axis per call. Recursion on the tuple rather than a -# loop over `1:D`, for the reason of `_solve_axes!`: the axes need not have the same -# quadrature type. -_scale_axes!(A, dims, ::Tuple{}, k) = A -function _scale_axes!(A, dims, qs::Tuple, k) - _scale_along!(A, dims, quadrature_weights(first(qs)), k) - _scale_axes!(A, dims, Base.tail(qs), k + 1) -end - +# The contraction, one axis per call. Recursion on the tuple rather than a loop over `1:D`, +# for the reason of `_solve_axes!`: the axes need not have the same quadrature type. # `src` holds an array of `n` entries whose first axis is the quadrature axis of `first(qs)`. # Each step writes `Φ * M` into `dst` transposed, as a (rest, N_k) array, which cycles the # contracted axis to the back; after the last step the axes are back in order, and the @@ -930,15 +930,22 @@ function l2_projection(q::TensorProductQuadrature{T, D}, f) where {T, D} end function l2_projection!(û::AbstractArray, q::TensorProductQuadrature{T, D}, f) where {T, D} + # The size of `û` is checked before `f` is sampled on the whole grid. + _check_coefficients(û, q) F = f isa AbstractArray ? f : quadrature_sample(q, f) l2_projection!(û, q, F, contract_workspace(q, eltype(û))) end function l2_projection!(û::AbstractArray, q::TensorProductQuadrature{T, D}, F::AbstractArray, work::ContractWorkspace) where {T, D} - size(û) == size(q.basis) || throw(DimensionMismatch( - "the coefficient array is $(size(û)) but the basis is $(size(q.basis))")) + _check_coefficients(û, q) contract!(û, q, F, work) mass_solve!(û, mass_operator(q), û, work.mass) return û end + +function _check_coefficients(û::AbstractArray, q::TensorProductQuadrature) + size(û) == size(q.basis) || throw(DimensionMismatch( + "the coefficient array is $(size(û)) but the basis is $(size(q.basis))")) + return nothing +end diff --git a/test/integration/workspaces.jl b/test/integration/workspaces.jl index 8abff93..4b3b7f6 100644 --- a/test/integration/workspaces.jl +++ b/test/integration/workspaces.jl @@ -245,6 +245,7 @@ end @test_throws ArgumentError mass_solve!(y, op, x, badsize) @test startswith(message(() -> mass_solve!(y, op, x, badsize)), "ArgumentError: the workspace has buffers of lengths " * + "$(map(n -> size(n, 1), mass_factors(other))) and " * "$(map(n -> size(n, 1), mass_factors(other))) but the operator needs " * "$(map(n -> size(n, 1), mass_factors(op)))") # a workspace of another dimension @@ -252,12 +253,24 @@ end b19 = BSplineBasis(UniformMesh(16, T(0) .. T(1)), 3) flat = mass_workspace(mass_operator(TensorProductQuadrature(b19 ⊗ b19 ⊗ b19 ⊗ b19)), T) @test_throws ArgumentError mass_solve!(y, op, x, flat) - # a long `dst` buffer - longdst = SimpleSplines.KroneckerMassWorkspace(work.src, - map(v -> [v; zero(T)], work.dst)) + # a long `src` buffer, a long `dst` buffer, and both too long by different amounts, + # which the message names each at its own place + dims = map(n -> size(n, 1), mass_factors(op)) + grow(vs, m) = map(v -> [v; zeros(T, m)], vs) + longsrc = SimpleSplines.KroneckerMassWorkspace(grow(work.src, 1), work.dst) + @test_throws ArgumentError mass_solve!(y, op, x, longsrc) + @test startswith(message(() -> mass_solve!(y, op, x, longsrc)), + "ArgumentError: the workspace has buffers of lengths $(dims .+ 1) and $(dims) " * + "but the operator needs $(dims)") + longdst = SimpleSplines.KroneckerMassWorkspace(work.src, grow(work.dst, 1)) @test_throws ArgumentError mass_solve!(y, op, x, longdst) @test startswith(message(() -> mass_solve!(y, op, x, longdst)), - "ArgumentError: the workspace has buffers of lengths") + "ArgumentError: the workspace has buffers of lengths $(dims) and $(dims .+ 1) " * + "but the operator needs $(dims)") + longboth = SimpleSplines.KroneckerMassWorkspace(grow(work.src, 2), grow(work.dst, 1)) + @test startswith(message(() -> mass_solve!(y, op, x, longboth)), + "ArgumentError: the workspace has buffers of lengths $(dims .+ 2) and " * + "$(dims .+ 1) but the operator needs $(dims)") # an array of the wrong length long = zeros(T, length(x) + 1) @test_throws DimensionMismatch mass_solve!(long, op, vec(x), work) @@ -318,6 +331,14 @@ end @test startswith(message(() -> l2_projection!(Lbad, q, F, work)), "DimensionMismatch: the coefficient array is $(size(Lbad)) but the basis is " * "$(size(B))") + # the three-argument form checks `û` before it samples `f` + calls = Ref(0) + counted = x -> (calls[] += 1; one(T)) + @test_throws DimensionMismatch l2_projection!(Lbad, q, counted) + @test calls[] == 0 + @test startswith(message(() -> l2_projection!(Lbad, q, counted)), + "DimensionMismatch: the coefficient array is $(size(Lbad)) but the basis is " * + "$(size(B))") # a workspace of the wrong element type names both types S = T === Float64 ? Float32 : Float64 @@ -331,13 +352,28 @@ end badsize = contract_workspace(TensorProductQuadrature(larger(B)), T) @test_throws ArgumentError contract!(L, q, F, badsize) @test startswith(message(() -> contract!(L, q, F, badsize)), - "ArgumentError: the workspace holds buffers of $(length(badsize.a)) entries but " * - "this quadrature needs $(length(work.a))") + "ArgumentError: the workspace holds buffers of $(length(badsize.a)) and " * + "$(length(badsize.b)) entries but this quadrature needs $(length(work.a))") @test length(badsize.a) != length(work.a) @test_throws ArgumentError l2_projection!(û, q, F, badsize) - # a long second buffer + # a long first buffer, a long second buffer, and both too long by different + # amounts, which the message names each at its own place + n = length(work.a) + longa = SimpleSplines.ContractWorkspace([work.a; zero(T)], work.b, work.mass) + @test_throws ArgumentError contract!(L, q, F, longa) + @test startswith(message(() -> contract!(L, q, F, longa)), + "ArgumentError: the workspace holds buffers of $(n + 1) and $(n) entries but " * + "this quadrature needs $(n)") longb = SimpleSplines.ContractWorkspace(work.a, [work.b; zero(T)], work.mass) @test_throws ArgumentError contract!(L, q, F, longb) + @test startswith(message(() -> contract!(L, q, F, longb)), + "ArgumentError: the workspace holds buffers of $(n) and $(n + 1) entries but " * + "this quadrature needs $(n)") + longab = SimpleSplines.ContractWorkspace([work.a; zeros(T, 2)], [work.b; zero(T)], + work.mass) + @test startswith(message(() -> contract!(L, q, F, longab)), + "ArgumentError: the workspace holds buffers of $(n + 2) and $(n + 1) entries " * + "but this quadrature needs $(n)") # a mass workspace of a larger operator wrongmass = SimpleSplines.ContractWorkspace(work.a, work.b, badsize.mass) @test_throws ArgumentError l2_projection!(û, q, F, wrongmass) From 1da79aae3471aba82d8fcd5e21da8918aa51656f Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 03:41:45 +0100 Subject: [PATCH 03/12] Describe the index arithmetic of _splitdims, which no longer feeds a reshape The per-axis loops read the array through linear indices, so the comment names those, and drops a measurement of the removed reshape form. Co-Authored-By: Claude Opus 5.5 --- src/tensorproduct.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 42734ca..5bcba2f 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -402,15 +402,14 @@ The tuple of one-dimensional [`MassOperator`](@ref)s, axis order. """ mass_factors(op::KroneckerMass) = op.ops -# The extents of a D-dimensional array before and after axis k, for the three-way reshape the -# per-axis loops below work on. +# The extents of a D-dimensional array before and after axis k, for the three-way +# (before, N_k, after) reading of the array that the per-axis loops below work on. # # Written as a loop rather than as `prod(dims[1:(k-1)]; init = 1)`. Slicing a tuple with a # runtime `k` gives a tuple whose *length* is not known to inference, so the product infers as -# `Any`, the `reshape` built from it infers as `Any`, and every element access in the loops -# that follow becomes a dynamic dispatch. Indexing the tuple one element at a time keeps `k` -# a run-time value and the result an `Int`. Measured on a 160x160 array, the difference in -# `_scale_along!` is 1862 us and 2.9 MB against 7.5 us and no allocation. +# `Any`, the linear indices built from it infer as `Any`, and every element access in the +# loops that follow becomes a dynamic dispatch. Indexing the tuple one element at a time keeps +# `k` a run-time value and the result an `Int`. @inline function _splitdims(dims::NTuple{D, Int}, k::Integer) where {D} before = 1 after = 1 From 1a8bad43f6c4c07fed225179331f51f388130aae Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 03:58:43 +0100 Subject: [PATCH 04/12] Correct the text of the contract! docstring, K8 and K9, and a stale reshape comment MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The contract! docstring says that only the size errors for û and F are those of the allocating forms; K8 names the four-argument mass_solve! that holds the copy; K9 gives the signature of evaluate_all! as it is; the _apply_along! comment no longer names a reshape; a CHANGELOG paragraph is rewrapped. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 10 +++++----- KNOWN_ISSUES.md | 6 ++++-- src/tensorproduct.jl | 10 +++++----- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6f79243..81de07d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -34,11 +34,11 @@ cells. Each pair is `Float64` / `Float32`. The solve along the axes of a `KroneckerMass` and the contraction along the axes of `contract` are unrolled over the tuple of axes, so an operator or a quadrature whose axes have different types no longer dispatches dynamically: that is the drop of the three-argument `mass_solve!` on -the periodic × clamped basis. The allocating `contract` and `l2_projection!` -allocate more than before, because each builds two buffers of the largest intermediate array, -and `l2_projection!` also the buffers of the mass solve. `mass_solve!` on a `KroneckerMass` -checks the length of `y` before it copies `x`, and raises a `DimensionMismatch` for a `y` of the -wrong length. Where the element types of the arguments are equal, the results of the allocating +the periodic × clamped basis. The allocating `contract` and `l2_projection!` allocate more +than before, because each builds two buffers of the largest intermediate array, and +`l2_projection!` also the buffers of the mass solve. `mass_solve!` on a `KroneckerMass` checks +the length of `y` before it copies `x`, and raises a `DimensionMismatch` for a `y` of the wrong +length. Where the element types of the arguments are equal, the results of the allocating forms are bit-for-bit those of `cb82e6f`. The three-argument `l2_projection!` now contracts in `eltype(û)`, not in the promoted type of the quadrature and the sample, so where these types differ its result can differ from that of `cb82e6f` by the round-off of the narrower type. diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index df27dc7..32e6823 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -101,7 +101,8 @@ ### K8 · `mass_solve!` on a `KroneckerMass` does not check the length of `x` -- **location:** `src/tensorproduct.jl`, `mass_solve!(y, op::KroneckerMass, x)` +- **location:** `src/tensorproduct.jl`, `mass_solve!(y, op::KroneckerMass, x, work)`, which + the three-argument form calls - **evidence:** the solve copies `x` into `y` with `copyto!`, which accepts a shorter `x`, and then solves on the whole of `y`. On the basis of K7, with `y = zeros(size(B)...)` (88 entries) and `x = ones(83)`, `mass_solve!(y, mass_operator(q), x)` returns with no error, on `cb82e6f` @@ -115,6 +116,7 @@ - **location:** `src/tensorproduct.jl`, `evaluate!(bufs, B::TensorProductBasis, û, x, d)` - **evidence:** the method takes `bufs` as an `NTuple{D, <:AbstractVector{R}}`, so a tuple of a `Vector` and a `SubArray`, or of vectors of different element types, raises a - `MethodError`. `evaluate_all!(bufs, B, x, d)` on `cb82e6f` has the same signature. + `MethodError`. `evaluate_all!(bufs, B, x, d)`, as on `cb82e6f`, takes `bufs` as an + `NTuple{D, <:AbstractVector}` and refuses the same tuples. - **kind:** found late - **found:** 2026-11-02 diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 5bcba2f..321b990 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -426,7 +426,7 @@ end # than taken from a `reshape`, because a `reshape` of an `Array` allocates its header. # # The fibre is copied into a contiguous buffer rather than passed as a view. A view into the -# middle index of a three-way reshape has stride `before`, and an FFTW plan or a LAPACK banded +# middle index of a three-way reading has stride `before`, and an FFTW plan or a LAPACK banded # solve addresses only contiguous memory -- so the one-dimensional solve would stage such a # fibre through its own buffer in any case. The two buffers `src` and `dst` have N_k entries # and are reused across the fibres of the axis. They can come from a caller's workspace, so @@ -789,10 +789,10 @@ The forms of [`contract`](@ref) and [`l2_projection!`](@ref) that take a workspa `l2_projection!` writes `l2_projection(q, F)` into `û`. `F` is the sample on the quadrature grid and is not changed. The workspace holds two buffers for the intermediate arrays of the contraction and a [`mass_workspace`](@ref) of `mass_operator(q)` for the solve. Wrong sizes -of `L`, `û` or `F` raise the `DimensionMismatch` of the allocating forms; a workspace whose -element type is not `eltype(L)`, or `eltype(û)`, or whose sizes do not fit `q`, raises an -`ArgumentError`. `contract` builds the two buffers on each call, and the three-argument -`l2_projection!` a whole workspace. +of `L`, `û` or `F` raise a `DimensionMismatch`, for `û` and `F` that of the allocating +forms; a workspace whose element type is not `eltype(L)`, or `eltype(û)`, or whose sizes do +not fit `q`, raises an `ArgumentError`. `contract` builds the two buffers on each call, and +the three-argument `l2_projection!` a whole workspace. A workspace belongs to one task. It does not make the quadrature safe for two tasks: a [`CirculantMass`](@ref) factor of its mass operator holds a scratch vector of its own. From d5a623ffb82ed286f83fed444bce519f89abd535 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 04:04:09 +0100 Subject: [PATCH 05/12] Give the evaluate! header the buffer type its method takes, and drop a history word Co-Authored-By: Claude Opus 5.5 --- src/tensorproduct.jl | 2 +- test/integration/workspaces.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 321b990..a08d1db 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -218,7 +218,7 @@ function evaluate(B::TensorProductBasis{T, D}, û::AbstractArray{S, D}, x, end @doc raw""" - evaluate!(bufs::NTuple{D,AbstractVector}, B::TensorProductBasis, û::AbstractArray, x, + evaluate!(bufs::NTuple{D,<:AbstractVector}, B::TensorProductBasis, û::AbstractArray, x, d = ntuple(_ -> 0, D)) The value of [`evaluate`](@ref)`(B, û, x, d)`, computed in the caller's per-axis buffers, diff --git a/test/integration/workspaces.jl b/test/integration/workspaces.jl index 4b3b7f6..94ac956 100644 --- a/test/integration/workspaces.jl +++ b/test/integration/workspaces.jl @@ -122,7 +122,7 @@ end "basis, got $(n)") @test_throws DimensionMismatch evaluate!(zeros(n), Spline(b, û), 0.37) end - # the coefficient check of `evaluate`, unchanged + # the coefficient check of `evaluate` @test_throws DimensionMismatch evaluate!(zeros(w), b, zeros(nbasis(b) + 1), 0.37) @test startswith(message(() -> evaluate!(zeros(w), b, zeros(nbasis(b) + 1), 0.37)), "DimensionMismatch: the coefficient vector has $(nbasis(b) + 1) entries") From 9f51b39330f70d0a111609a8a84ca60b7e2fd5ad Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 09:59:14 +0100 Subject: [PATCH 06/12] Bind the type parameters of evaluate! and KroneckerMassWorkspace, check x in mass_solve!, refuse narrow buffers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - evaluate! on a TensorProductBasis takes NTuple{D, <:AbstractVector} and reads the element type from the buffers; KroneckerMassWorkspace keeps only its parametric constructor. Aqua's unbound-argument check passes on Julia 1.14. - mass_solve!(y, op::KroneckerMass, x, work) raises a DimensionMismatch for an x whose length is not that of y (K8, also in 0.3.4). - evaluate! raises an ArgumentError for a buffer that cannot hold the value, rather than rounding to it. - The four-argument l2_projection! takes only a workspace with a mass workspace, so it no longer writes û before a MethodError. - _check_coefficients takes the quadrature first, as its siblings do. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 31 +++++++++++++++++++---------- KNOWN_ISSUES.md | 18 +++-------------- src/bspline.jl | 11 +++++++---- src/tensorproduct.jl | 26 ++++++++++++++++-------- test/integration/workspaces.jl | 36 +++++++++++++++++++++++++++++++--- 5 files changed, 82 insertions(+), 40 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 81de07d..c9c5a99 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,8 +15,10 @@ on a one-dimensional basis with a buffer of `local_width(b)` entries; `evaluate! d)` does it on a `TensorProductBasis` with one buffer per axis. `mass_workspace(op, T)` and `mass_solve!(y, op, x, work)` solve with a `KroneckerMass`; `contract_workspace(q, T)`, `contract!(L, q, F, work, d)` and `l2_projection!(û, q, F, work)` contract and project on a -`TensorProductQuadrature`. Each buffer and workspace belongs to one task. The allocating forms -stay, and build their buffers and call these. +`TensorProductQuadrature`. Each buffer and workspace belongs to one task. A buffer whose element +type cannot hold the promoted type of the basis, the coefficients and `x` raises an +`ArgumentError`, and is not rounded to. The allocating forms stay, and build their buffers and +call these. Bytes per call, measured with `@allocated` after a warm call through a function barrier, in a cold process on Julia 1.13.3, on `cb82e6f` and on this release. The one-dimensional basis is @@ -34,18 +36,27 @@ cells. Each pair is `Float64` / `Float32`. The solve along the axes of a `KroneckerMass` and the contraction along the axes of `contract` are unrolled over the tuple of axes, so an operator or a quadrature whose axes have different types no longer dispatches dynamically: that is the drop of the three-argument `mass_solve!` on -the periodic × clamped basis. The allocating `contract` and `l2_projection!` allocate more -than before, because each builds two buffers of the largest intermediate array, and -`l2_projection!` also the buffers of the mass solve. `mass_solve!` on a `KroneckerMass` checks -the length of `y` before it copies `x`, and raises a `DimensionMismatch` for a `y` of the wrong -length. Where the element types of the arguments are equal, the results of the allocating -forms are bit-for-bit those of `cb82e6f`. The three-argument `l2_projection!` now contracts in -`eltype(û)`, not in the promoted type of the quadrature and the sample, so where these types -differ its result can differ from that of `cb82e6f` by the round-off of the narrower type. +the periodic × clamped basis. The allocating `contract` and `l2_projection!` allocate more than +before, because each builds two buffers of the largest intermediate array, and `l2_projection!` +also the buffers of the mass solve. `mass_solve!` on a `KroneckerMass` checks the lengths of `y` +and `x` before it copies `x`, and raises a `DimensionMismatch` for a `y` or an `x` whose length +is not that of the operator. Where the element types of the arguments are equal, the results of +the allocating forms are bit-for-bit those of `cb82e6f`. The three-argument `l2_projection!` now +contracts in `eltype(û)`, not in the promoted type of the quadrature and the sample, so where +these types differ its result can differ from that of `cb82e6f` by the round-off of the narrower +type. `ForwardDiff` joins the test dependencies: the tests evaluate with `ForwardDiff.Dual` arguments and buffers, which also allocate nothing. +### Bug Fixes + +`mass_solve!` on a `KroneckerMass`, in the three-argument form and so in `ldiv!` and `\`, accepted +an `x` shorter than `y` in 0.3.4, with no error. It copied `x` into the start of `y` and solved the +whole of `y`, so the solve ran on whatever the rest of `y` held. It now raises a +`DimensionMismatch` for an `x` whose length is not that of `y`. This fixes known issue K8. Earlier +releases were not checked. + ## [0.3.4] — 2026-10-31 ### Bug Fixes diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 32e6823..c4b70b1 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -99,24 +99,12 @@ - **kind:** defect - **found:** 2026-11-02; older than `cb82e6f` -### K8 · `mass_solve!` on a `KroneckerMass` does not check the length of `x` - -- **location:** `src/tensorproduct.jl`, `mass_solve!(y, op::KroneckerMass, x, work)`, which - the three-argument form calls -- **evidence:** the solve copies `x` into `y` with `copyto!`, which accepts a shorter `x`, and - then solves on the whole of `y`. On the basis of K7, with `y = zeros(size(B)...)` (88 entries) - and `x = ones(83)`, `mass_solve!(y, mass_operator(q), x)` returns with no error, on `cb82e6f` - and on the branch that adds caller-held workspaces alike; the last five entries of `x` are - taken from the zeros of `y`. -- **kind:** defect -- **found:** 2026-11-02; older than `cb82e6f` - ### K9 · `evaluate!` on a `TensorProductBasis` takes buffers of one type only - **location:** `src/tensorproduct.jl`, `evaluate!(bufs, B::TensorProductBasis, û, x, d)` -- **evidence:** the method takes `bufs` as an `NTuple{D, <:AbstractVector{R}}`, so a tuple +- **evidence:** the method takes `bufs` as an `NTuple{D, <:AbstractVector}`, so a tuple of a `Vector` and a `SubArray`, or of vectors of different element types, raises a - `MethodError`. `evaluate_all!(bufs, B, x, d)`, as on `cb82e6f`, takes `bufs` as an - `NTuple{D, <:AbstractVector}` and refuses the same tuples. + `MethodError`. `evaluate_all!(bufs, B, x, d)`, which it calls, takes `bufs` as the same + type, as on `cb82e6f`, and refuses the same tuples. - **kind:** found late - **found:** 2026-11-02 diff --git a/src/bspline.jl b/src/bspline.jl index 373226d..47683b1 100644 --- a/src/bspline.jl +++ b/src/bspline.jl @@ -639,8 +639,8 @@ caller's buffer `buf`, which allocates nothing. `buf` must have [`local_width`](@ref)`(b)` entries, or a `DimensionMismatch` is raised; its contents on entry do not matter, and on return they are the local block of [`evaluate_all!`](@ref). Its element type is the type the sum is accumulated in, so it must -hold `promote_type(eltype(b), eltype(û), typeof(x))`: a `ForwardDiff.Dual` argument takes a -buffer of `Dual`s. +hold `promote_type(eltype(b), eltype(û), typeof(x))`, or an `ArgumentError` is raised: a +`ForwardDiff.Dual` argument takes a buffer of `Dual`s. A buffer belongs to one task. Two tasks that evaluate at once need a buffer each. @@ -653,14 +653,17 @@ julia> evaluate!(buf, b, û, 0.3, 1) == evaluate(b, û, 0.3, 1) true ``` """ -function evaluate!(buf::AbstractVector, b::AbstractBSplineBasis, û::AbstractVector, - x::Number, d::Integer = 0) +function evaluate!(buf::AbstractVector{B}, b::AbstractBSplineBasis{T}, + û::AbstractVector{S}, x::Number, d::Integer = 0) where {B, T, S} length(û) == nbasis(b) || throw(DimensionMismatch( "the coefficient vector has $(length(û)) entries but the basis has $(nbasis(b))")) w = local_width(b) length(buf) == w || throw(DimensionMismatch( "evaluate! needs a buffer of $(w) entries for this basis, got $(length(buf)); " * "see `local_width`")) + R = _evaltype(promote_type(T, S), typeof(x)) + promote_type(B, R) === B || throw(ArgumentError( + "evaluate! needs a buffer that holds $(R), got a buffer of $(B)")) _evaluate_block!(buf, b, û, x, d) end diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index a08d1db..6c45df6 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -227,17 +227,21 @@ which allocates nothing. `bufs[k]` must have `local_width(bases(B)[k])` entries, or a `DimensionMismatch` is raised; their contents on entry do not matter, and on return they are the blocks of [`evaluate_all!`](@ref). Their element type is the type the sum is accumulated in, so it must -hold the promoted type of the basis, the coefficients and `x`. +hold the promoted type of the basis, the coefficients and `x`, or an `ArgumentError` is raised. The buffers belong to one task. Two tasks that evaluate at once need a tuple each. """ -function evaluate!(bufs::NTuple{D, <:AbstractVector{R}}, B::TensorProductBasis{T, D}, - û::AbstractArray{S, D}, x, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {R, T, D, S} +function evaluate!(bufs::NTuple{D, <:AbstractVector}, B::TensorProductBasis{T, D}, + û::AbstractArray{S, D}, x, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {T, D, S} size(û) == size(B) || throw(DimensionMismatch( "the coefficient array is $(size(û)) but the basis is $(size(B))")) widths = local_width(B) map(length, bufs) == widths || throw(DimensionMismatch( "evaluate! needs buffers of $(widths) entries on the axes, got $(map(length, bufs))")) + R = eltype(first(bufs)) + Rₓ = _evaltype(promote_type(T, S), eltype(x)) + promote_type(R, Rₓ) === R || throw(ArgumentError( + "evaluate! needs buffers that hold $(Rₓ), got buffers of $(R)")) # Only the local block contributes, so this is the fast path rather than a sum over the # whole array. `evaluate_all!` gives the p_d+1 nonzero factors on each axis; their outer @@ -471,6 +475,9 @@ end struct KroneckerMassWorkspace{T, D} src::NTuple{D, Vector{T}} dst::NTuple{D, Vector{T}} + + # Only with its parameters: at `D = 0` the arguments would leave `T` unbound. + KroneckerMassWorkspace{T, D}(src, dst) where {T, D} = new{T, D}(src, dst) end Base.eltype(::KroneckerMassWorkspace{T}) where {T} = T @@ -484,7 +491,7 @@ which allocates nothing. The workspace holds two buffers of `size(mass_factors(op)[k], 1)` entries for each axis `k`. The solve takes `y` and `x` as arrays of the size of the tensor-product basis, or as vectors of -`size(op, 1)` entries, and `y === x` is allowed. A `y` of another length raises a +`size(op, 1)` entries, and `y === x` is allowed. A `y` or an `x` of another length raises a `DimensionMismatch`; a workspace whose element type is not `eltype(y)`, or whose buffer lengths do not fit `op`, raises an `ArgumentError`. `mass_solve!(y, op, x)` builds a workspace of `eltype(y)` on each call. @@ -519,6 +526,8 @@ function mass_solve!(y::AbstractArray, op::KroneckerMass{T, D}, x::AbstractArray "$(map(length, work.dst)) but the operator needs $(op.dims)")) length(y) == prod(op.dims) || throw(DimensionMismatch( "the array has $(length(y)) entries but the operator acts on $(op.dims)")) + length(x) == length(y) || throw(DimensionMismatch( + "the argument has $(length(x)) entries but the operator acts on $(op.dims)")) y === x || copyto!(y, x) _solve_axes!(y, op.dims, op.ops, work.src, work.dst, 1) return y @@ -930,20 +939,21 @@ end function l2_projection!(û::AbstractArray, q::TensorProductQuadrature{T, D}, f) where {T, D} # The size of `û` is checked before `f` is sampled on the whole grid. - _check_coefficients(û, q) + _check_coefficients(q, û) F = f isa AbstractArray ? f : quadrature_sample(q, f) l2_projection!(û, q, F, contract_workspace(q, eltype(û))) end function l2_projection!(û::AbstractArray, q::TensorProductQuadrature{T, D}, - F::AbstractArray, work::ContractWorkspace) where {T, D} - _check_coefficients(û, q) + F::AbstractArray, + work::ContractWorkspace{<:Any, <:KroneckerMassWorkspace}) where {T, D} + _check_coefficients(q, û) contract!(û, q, F, work) mass_solve!(û, mass_operator(q), û, work.mass) return û end -function _check_coefficients(û::AbstractArray, q::TensorProductQuadrature) +function _check_coefficients(q::TensorProductQuadrature, û::AbstractArray) size(û) == size(q.basis) || throw(DimensionMismatch( "the coefficient array is $(size(û)) but the basis is $(size(q.basis))")) return nothing diff --git a/test/integration/workspaces.jl b/test/integration/workspaces.jl index 94ac956..663e026 100644 --- a/test/integration/workspaces.jl +++ b/test/integration/workspaces.jl @@ -126,6 +126,13 @@ end @test_throws DimensionMismatch evaluate!(zeros(w), b, zeros(nbasis(b) + 1), 0.37) @test startswith(message(() -> evaluate!(zeros(w), b, zeros(nbasis(b) + 1), 0.37)), "DimensionMismatch: the coefficient vector has $(nbasis(b) + 1) entries") + # a buffer that cannot hold the value, of a narrower type or for a Dual argument + @test_throws ArgumentError evaluate!(zeros(Float32, w), b, û, 0.37) + @test startswith(message(() -> evaluate!(zeros(Float32, w), b, û, 0.37)), + "ArgumentError: evaluate! needs a buffer that holds Float64, got a buffer of " * + "Float32") + @test_throws ArgumentError evaluate!(zeros(Float32, w), Spline(b, û), 0.37) + @test_throws ArgumentError evaluate!(zeros(w), b, û, dual(0.37)) end end @@ -207,6 +214,15 @@ end @test startswith(message(() -> evaluate!(bufs, B, ûbad, first(points))), "DimensionMismatch: the coefficient array is $(size(ûbad)) but the basis is " * "$(size(B))") + # buffers that cannot hold the value, for Dual arguments or of a narrower type + @test_throws ArgumentError evaluate!(bufs, B, û, map(dual, first(points))) + if T === Float64 + narrow = map(w -> zeros(Float32, w), local_width(B)) + @test_throws ArgumentError evaluate!(narrow, B, û, first(points)) + @test startswith(message(() -> evaluate!(narrow, B, û, first(points))), + "ArgumentError: evaluate! needs buffers that hold Float64, got buffers of " * + "Float32") + end end end @@ -257,17 +273,17 @@ end # which the message names each at its own place dims = map(n -> size(n, 1), mass_factors(op)) grow(vs, m) = map(v -> [v; zeros(T, m)], vs) - longsrc = SimpleSplines.KroneckerMassWorkspace(grow(work.src, 1), work.dst) + longsrc = typeof(work)(grow(work.src, 1), work.dst) @test_throws ArgumentError mass_solve!(y, op, x, longsrc) @test startswith(message(() -> mass_solve!(y, op, x, longsrc)), "ArgumentError: the workspace has buffers of lengths $(dims .+ 1) and $(dims) " * "but the operator needs $(dims)") - longdst = SimpleSplines.KroneckerMassWorkspace(work.src, grow(work.dst, 1)) + longdst = typeof(work)(work.src, grow(work.dst, 1)) @test_throws ArgumentError mass_solve!(y, op, x, longdst) @test startswith(message(() -> mass_solve!(y, op, x, longdst)), "ArgumentError: the workspace has buffers of lengths $(dims) and $(dims .+ 1) " * "but the operator needs $(dims)") - longboth = SimpleSplines.KroneckerMassWorkspace(grow(work.src, 2), grow(work.dst, 1)) + longboth = typeof(work)(grow(work.src, 2), grow(work.dst, 1)) @test startswith(message(() -> mass_solve!(y, op, x, longboth)), "ArgumentError: the workspace has buffers of lengths $(dims .+ 2) and " * "$(dims .+ 1) but the operator needs $(dims)") @@ -277,6 +293,15 @@ end @test startswith(message(() -> mass_solve!(long, op, vec(x), work)), "DimensionMismatch: the array has $(length(x) + 1) entries but the operator " * "acts on $(size(x))") + # an argument shorter or longer than `y`, through both forms + for n in (length(x) - 5, length(x) + 1) + xbad = ones(T, n) + @test_throws DimensionMismatch mass_solve!(y, op, xbad, work) + @test_throws DimensionMismatch mass_solve!(y, op, xbad) + @test startswith(message(() -> mass_solve!(y, op, xbad, work)), + "DimensionMismatch: the argument has $(n) entries but the operator acts on " * + "$(size(x))") + end end end @@ -379,5 +404,10 @@ end @test_throws ArgumentError l2_projection!(û, q, F, wrongmass) @test startswith(message(() -> l2_projection!(û, q, F, wrongmass)), "ArgumentError: the workspace has buffers of lengths") + # a workspace without a mass workspace, as `contract` builds it, leaves `û` alone + nomass = SimpleSplines.ContractWorkspace(work.a, work.b, nothing) + û₀ = copy(û) + @test_throws MethodError l2_projection!(û, q, F, nomass) + @test û == û₀ end end From 9e2f5c86fba00eeaa7db85a438d2ebfd3dc98249 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 11:37:44 +0100 Subject: [PATCH 07/12] Qualify the allocation-free claim and the task-sharing rule A periodic axis on a non-uniform mesh solves with a FactorizedMass, which allocates a CHOLMOD temporary. A shared mass operator or quadrature is unsafe across tasks because a CirculantMass factor writes scratch space. Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 21 ++++++++++++++------- src/tensorproduct.jl | 9 ++++++--- 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c9c5a99..abf43a3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,15 +10,22 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### New Features The hot evaluation, solve and projection paths take buffers that the caller holds, and then -allocate nothing. `evaluate!(buf, b, û, x, d)` and `evaluate!(buf, s, x, d)` evaluate a spline -on a one-dimensional basis with a buffer of `local_width(b)` entries; `evaluate!(bufs, B, û, x, -d)` does it on a `TensorProductBasis` with one buffer per axis. `mass_workspace(op, T)` and +allocate nothing, except for the solve on a periodic axis with a non-uniform mesh +(`GradedMesh` or `RandomMesh`). The `FactorizedMass` of that axis allocates a CHOLMOD +temporary, because CHOLMOD has no in-place `ldiv!`. Measured in a cold process with `Float64`, +periodic × clamped, cubic, on 8 × 8 cells: 0 bytes on a `UniformMesh`, and 9328 bytes per +`mass_solve!(y, op, x, work)` and per `l2_projection!(û, q, F, work)` on a `GradedMesh` or a +`RandomMesh`. +`evaluate!(buf, b, û, x, d)` and `evaluate!(buf, s, x, d)` evaluate a spline on a +one-dimensional basis with a buffer of `local_width(b)` entries; `evaluate!(bufs, B, û, x, d)` +does it on a `TensorProductBasis` with one buffer per axis. `mass_workspace(op, T)` and `mass_solve!(y, op, x, work)` solve with a `KroneckerMass`; `contract_workspace(q, T)`, `contract!(L, q, F, work, d)` and `l2_projection!(û, q, F, work)` contract and project on a -`TensorProductQuadrature`. Each buffer and workspace belongs to one task. A buffer whose element -type cannot hold the promoted type of the basis, the coefficients and `x` raises an -`ArgumentError`, and is not rounded to. The allocating forms stay, and build their buffers and -call these. +`TensorProductQuadrature`. Each buffer and workspace belongs to one task. Tasks must not share a +mass operator or `TensorProductQuadrature`, because a `CirculantMass` factor writes scratch space +during each solve. A buffer whose element type cannot hold the promoted type of the basis, the +coefficients and `x` raises an `ArgumentError`, and is not rounded to. The allocating forms stay, +and build their buffers and call these. Bytes per call, measured with `@allocated` after a warm call through a function barrier, in a cold process on Julia 1.13.3, on `cb82e6f` and on this release. The one-dimensional basis is diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 6c45df6..7f63e18 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -487,7 +487,8 @@ Base.eltype(::KroneckerMassWorkspace{T}) where {T} = T mass_solve!(y, op::KroneckerMass, x, work) A workspace of element type `T` for [`mass_solve!`](@ref) on `op`, and the solve that uses it, -which allocates nothing. +which allocates nothing unless an axis of `op` is a [`FactorizedMass`](@ref). The solve of that +axis, a periodic basis on a non-uniform mesh, allocates a CHOLMOD temporary. The workspace holds two buffers of `size(mass_factors(op)[k], 1)` entries for each axis `k`. The solve takes `y` and `x` as arrays of the size of the tensor-product basis, or as vectors of @@ -792,7 +793,8 @@ end l2_projection!(û, q::TensorProductQuadrature, F::AbstractArray, work) The forms of [`contract`](@ref) and [`l2_projection!`](@ref) that take a workspace from -[`contract_workspace`](@ref), and allocate nothing. +[`contract_workspace`](@ref), and allocate nothing. The exception is a periodic axis on a +non-uniform mesh: the solve of its [`FactorizedMass`](@ref) allocates a CHOLMOD temporary. `contract!` writes [`contract`](@ref)`(q, F, d)` into `L`, of size `size(basis(q))`; `l2_projection!` writes `l2_projection(q, F)` into `û`. `F` is the sample on the quadrature @@ -910,7 +912,8 @@ The coefficient array of the ``L^2`` projection of `f` onto the tensor-product s `f` may be a function of a `D`-tuple of coordinates, or an array already sampled on the quadrature grid — see [`quadrature_sample`](@ref). The result has size `size(basis(q))`. The form with a workspace from [`contract_workspace`](@ref) takes the sample and allocates -nothing. +nothing, except on a periodic axis with a non-uniform mesh, where the solve of its +[`FactorizedMass`](@ref) allocates a CHOLMOD temporary. ```jldoctest julia> B = BSplineBasis(UniformMesh(16, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(16, 0 .. 1), 3); From 93c41e7342baf75c4d98871cec9eaf496d1a46cd Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 11:55:49 +0100 Subject: [PATCH 08/12] evaluate!: the buffer type follows the evaluation point A Float64 basis with Float64 coefficients at a Float32 point now works in a Float32 buffer, returns a Float32 and allocates nothing. The buffer must hold the type of the point; a Float64 point into a Float32 buffer and a Dual point into a Float64 buffer still raise an ArgumentError. The triangular scheme of evaluate_all! runs in the real type of the buffer. Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 9 +++-- src/bspline.jl | 31 ++++++++++------ src/tensorproduct.jl | 13 +++---- test/integration/workspaces.jl | 64 ++++++++++++++++++++++++++++++---- 4 files changed, 92 insertions(+), 25 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index abf43a3..c3fd5f7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,8 +23,13 @@ does it on a `TensorProductBasis` with one buffer per axis. `mass_workspace(op, `contract!(L, q, F, work, d)` and `l2_projection!(û, q, F, work)` contract and project on a `TensorProductQuadrature`. Each buffer and workspace belongs to one task. Tasks must not share a mass operator or `TensorProductQuadrature`, because a `CirculantMass` factor writes scratch space -during each solve. A buffer whose element type cannot hold the promoted type of the basis, the -coefficients and `x` raises an `ArgumentError`, and is not rounded to. The allocating forms stay, +during each solve. The element type of the buffer follows the evaluation point, not the basis or the coefficients. +`evaluate!` computes in the buffer's type and returns it. The buffer must hold `typeof(x)` +(`eltype(x)` for the tensor-product form), or an `ArgumentError` is raised. So a `Float64` basis +with `Float64` coefficients at a `Float32` point works in a `Float32` buffer. It returns a +`Float32` and allocates 0 bytes, measured with `@allocated` after a warm call on a cubic clamped +basis. A `Float64` point into a `Float32` buffer raises, and so does a `ForwardDiff.Dual` point +into a `Float64` buffer. The allocating forms stay, and build their buffers and call these. Bytes per call, measured with `@allocated` after a warm call through a function barrier, in a diff --git a/src/bspline.jl b/src/bspline.jl index 47683b1..c38b087 100644 --- a/src/bspline.jl +++ b/src/bspline.jl @@ -638,9 +638,11 @@ caller's buffer `buf`, which allocates nothing. `buf` must have [`local_width`](@ref)`(b)` entries, or a `DimensionMismatch` is raised; its contents on entry do not matter, and on return they are the local block of -[`evaluate_all!`](@ref). Its element type is the type the sum is accumulated in, so it must -hold `promote_type(eltype(b), eltype(û), typeof(x))`, or an `ArgumentError` is raised: a -`ForwardDiff.Dual` argument takes a buffer of `Dual`s. +[`evaluate_all!`](@ref). Its element type is the type the sum is computed in, and the type of +the result. It follows the point, not the basis or the coefficients: it must hold `typeof(x)`, +or an `ArgumentError` is raised. A `Float64` basis with `Float64` coefficients at a `Float32` +point takes a `Float32` buffer and returns a `Float32`. A `ForwardDiff.Dual` point takes a +buffer of `Dual`s, and a `Float64` point does not take a `Float32` buffer. A buffer belongs to one task. Two tasks that evaluate at once need a buffer each. @@ -661,9 +663,8 @@ function evaluate!(buf::AbstractVector{B}, b::AbstractBSplineBasis{T}, length(buf) == w || throw(DimensionMismatch( "evaluate! needs a buffer of $(w) entries for this basis, got $(length(buf)); " * "see `local_width`")) - R = _evaltype(promote_type(T, S), typeof(x)) - promote_type(B, R) === B || throw(ArgumentError( - "evaluate! needs a buffer that holds $(R), got a buffer of $(B)")) + promote_type(B, typeof(x)) === B || throw(ArgumentError( + "evaluate! needs a buffer that holds the point type $(typeof(x)), got a buffer of $(B)")) _evaluate_block!(buf, b, û, x, d) end @@ -688,7 +689,7 @@ function _evaluate_block!(values::AbstractVector{R}, b::AbstractBSplineBasis, # index outside 1:N means the block of this cell reaches past the end of the basis. # Such a term is genuinely absent rather than zero-valued, so it is skipped. 1 ≤ j ≤ nbasis(b) || continue - v += values[t] * û[j] + v += convert(R, values[t] * û[j]) end return v end @@ -870,6 +871,14 @@ function evaluate_all!(values::AbstractVector{R}, b::AbstractBSplineBasis, x::Nu q = p - Int(d) + # The arithmetic runs in the real type of the element type of `values`: the cell is found + # above in the type of the basis, and the knots and the point are converted once the cell + # is known. A `Float64` basis into a `Float32` buffer therefore does a `Float32` sum, with + # no `Float64` intermediate to convert back. `real`, so that a complex buffer keeps real + # knots and the comparisons of the derivative step stay real. + Q = real(R) + x̃ = convert(Q, x̃) + # De Boor's triangular scheme for the q+1 nonzero splines of degree q on span k. # values[1:q+1] holds φ_{k-q}^q .. φ_k^q afterwards. values[1] = one(R) @@ -878,8 +887,8 @@ function evaluate_all!(values::AbstractVector{R}, b::AbstractBSplineBasis, x::Nu for r in 1:j # right = kv[k+r] - x, left = kv[k+1-(j-r+1)] - ... written out below to keep # the index arithmetic of the standard form visible. - right = kv[k + r] - x̃ - left = x̃ - kv[k + 1 - j + r - 1] + right = convert(Q, kv[k + r]) - x̃ + left = x̃ - convert(Q, kv[k + 1 - j + r - 1]) temp = values[r] / (right + left) values[r] = saved + right * temp saved = left * temp @@ -899,8 +908,8 @@ function evaluate_all!(values::AbstractVector{R}, b::AbstractBSplineBasis, x::Nu below = t ≥ 2 ? values[t - 1] : zero(R) above = t ≤ count ? values[t] : zero(R) j = k - count - 1 + t - da = kv[j + r] - kv[j] - db = kv[j + r + 1] - kv[j + 1] + da = convert(Q, kv[j + r] - kv[j]) + db = convert(Q, kv[j + r + 1] - kv[j + 1]) out = zero(R) da > 0 && (out += below / da) db > 0 && (out -= above / db) diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 7f63e18..3da328d 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -226,8 +226,10 @@ which allocates nothing. `bufs[k]` must have `local_width(bases(B)[k])` entries, or a `DimensionMismatch` is raised; their contents on entry do not matter, and on return they are the blocks of -[`evaluate_all!`](@ref). Their element type is the type the sum is accumulated in, so it must -hold the promoted type of the basis, the coefficients and `x`, or an `ArgumentError` is raised. +[`evaluate_all!`](@ref). Their element type is the type the sum is computed in, and the type of +the result. It follows the point, not the basis or the coefficients: it must hold `eltype(x)`, +or an `ArgumentError` is raised. A `Float64` basis at a `Float32` point takes `Float32` +buffers and returns a `Float32`. The buffers belong to one task. Two tasks that evaluate at once need a tuple each. """ @@ -239,9 +241,8 @@ function evaluate!(bufs::NTuple{D, <:AbstractVector}, B::TensorProductBasis{T, D map(length, bufs) == widths || throw(DimensionMismatch( "evaluate! needs buffers of $(widths) entries on the axes, got $(map(length, bufs))")) R = eltype(first(bufs)) - Rₓ = _evaltype(promote_type(T, S), eltype(x)) - promote_type(R, Rₓ) === R || throw(ArgumentError( - "evaluate! needs buffers that hold $(Rₓ), got buffers of $(R)")) + promote_type(R, eltype(x)) === R || throw(ArgumentError( + "evaluate! needs buffers that hold the point type $(eltype(x)), got buffers of $(R)")) # Only the local block contributes, so this is the fast path rather than a sum over the # whole array. `evaluate_all!` gives the p_d+1 nonzero factors on each axis; their outer @@ -264,7 +265,7 @@ function evaluate!(bufs::NTuple{D, <:AbstractVector}, B::TensorProductBasis{T, D # The weight is formed here rather than accumulated inside the `ntuple` above: a # variable assigned from within that closure is boxed, which infers as `Any` and makes # every element access below a dynamic dispatch. - v += prod(ntuple(k -> bufs[k][t[k]], D)) * û[idx...] + v += convert(R, prod(ntuple(k -> bufs[k][t[k]], D)) * û[idx...]) end return v end diff --git a/test/integration/workspaces.jl b/test/integration/workspaces.jl index 663e026..8a2f2c0 100644 --- a/test/integration/workspaces.jl +++ b/test/integration/workspaces.jl @@ -126,13 +126,65 @@ end @test_throws DimensionMismatch evaluate!(zeros(w), b, zeros(nbasis(b) + 1), 0.37) @test startswith(message(() -> evaluate!(zeros(w), b, zeros(nbasis(b) + 1), 0.37)), "DimensionMismatch: the coefficient vector has $(nbasis(b) + 1) entries") - # a buffer that cannot hold the value, of a narrower type or for a Dual argument + # a buffer that cannot hold the type of the point, a Float64 point into a Float32 + # buffer or a Dual point into a Float64 buffer @test_throws ArgumentError evaluate!(zeros(Float32, w), b, û, 0.37) @test startswith(message(() -> evaluate!(zeros(Float32, w), b, û, 0.37)), - "ArgumentError: evaluate! needs a buffer that holds Float64, got a buffer of " * - "Float32") + "ArgumentError: evaluate! needs a buffer that holds the point type Float64, " * + "got a buffer of Float32") @test_throws ArgumentError evaluate!(zeros(Float32, w), Spline(b, û), 0.37) @test_throws ArgumentError evaluate!(zeros(w), b, û, dual(0.37)) + @test startswith(message(() -> evaluate!(zeros(w), b, û, dual(0.37))), + "ArgumentError: evaluate! needs a buffer that holds the point type $(Dual12), " * + "got a buffer of Float64") + end +end + +@testset "$(rpad("evaluate! of a Float64 spline at a Float32 point",76))" begin + for b in bases_1d(Float64) + û = randn(nbasis(b)) + s = Spline(b, û) + buf = zeros(Float32, local_width(b)) + for x in Float32[0.37, 0.5, 0, 1, prevfloat(1.0f0), -0.25, 1.75], d in 0:2 + + ref = Float32(evaluate(b, û, Float64(x), d)) + v = @inferred evaluate!(buf, b, û, x, d) + @test v isa Float32 + # The sum cancels, so the floor is the size of its terms: the coefficients, times + # the size of a basis function after d derivatives on 8 cells. + atol = 8 * eps(Float32) * maximum(abs, û) * (8 * degree(b))^d + @test isapprox(v, ref; rtol = 8 * eps(Float32), atol) + w = @inferred evaluate!(buf, s, x, d) + @test w isa Float32 + @test w == v + end + if MEASURE + @test allocs(evaluate!, buf, b, û, 0.37f0, 1) == 0 + @test allocs(evaluate!, buf, s, 0.37f0, 1) == 0 + end + end + + # Float64 point types are unaffected, and a wider buffer holds a Float32 point + b = bases_1d(Float64)[2] + û = randn(nbasis(b)) + @test evaluate!(zeros(local_width(b)), b, û, 0.37f0) isa Float64 + + for B in (basis_2d(Float64), basis_3d(Float64)) + D = ndims(B) + û = randn(size(B)...) + bufs = map(w -> zeros(Float32, w), local_width(B)) + x = D == 2 ? (0.37f0, 0.4f0) : (0.37f0, 1.2f0, -0.3f0) + for d in (ntuple(_ -> 0, D), ntuple(k -> k == 1 ? 1 : 0, D)) + ref = Float32(evaluate(B, û, Float64.(x), d)) + v = @inferred evaluate!(bufs, B, û, x, d) + @test v isa Float32 + atol = 8 * eps(Float32) * maximum(abs, û) * 24^sum(d) + @test isapprox(v, ref; rtol = 8 * eps(Float32), atol) + end + if MEASURE + d = ntuple(k -> k == 1 ? 1 : 0, D) + @test allocs(evaluate!, bufs, B, û, x, d) == 0 + end end end @@ -214,14 +266,14 @@ end @test startswith(message(() -> evaluate!(bufs, B, ûbad, first(points))), "DimensionMismatch: the coefficient array is $(size(ûbad)) but the basis is " * "$(size(B))") - # buffers that cannot hold the value, for Dual arguments or of a narrower type + # buffers that cannot hold the type of the point, for Dual points or of a narrower type @test_throws ArgumentError evaluate!(bufs, B, û, map(dual, first(points))) if T === Float64 narrow = map(w -> zeros(Float32, w), local_width(B)) @test_throws ArgumentError evaluate!(narrow, B, û, first(points)) @test startswith(message(() -> evaluate!(narrow, B, û, first(points))), - "ArgumentError: evaluate! needs buffers that hold Float64, got buffers of " * - "Float32") + "ArgumentError: evaluate! needs buffers that hold the point type Float64, " * + "got buffers of Float32") end end end From 38ec2f412a9165e8f1dd4c453f896fd10654e926 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 13:23:24 +0100 Subject: [PATCH 09/12] The allocating evaluate forms return the type of the point evaluate, evaluate_all and the Spline, tensor-product, recombined and polar forms take the result type from the point, as evaluate! does. Complex coefficients give a complex result in that type. A point that is not a floating-point type keeps the promotion of the basis and the point. Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 17 +++++- src/bspline.jl | 35 +++++++++--- src/polar.jl | 23 +++++--- src/recombined.jl | 2 +- src/tensorproduct.jl | 10 +++- test/integration/workspaces.jl | 99 ++++++++++++++++++++++++++++++++++ test/polar.jl | 46 ++++++++++++++++ 7 files changed, 216 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bcbf2de..d8ee1f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,7 +23,22 @@ does it on a `TensorProductBasis` with one buffer per axis. `mass_workspace(op, `contract!(L, q, F, work, d)` and `l2_projection!(û, q, F, work)` contract and project on a `TensorProductQuadrature`. Each buffer and workspace belongs to one task. Tasks must not share a mass operator or `TensorProductQuadrature`, because a `CirculantMass` factor writes scratch space -during each solve. The element type of the buffer follows the evaluation point, not the basis or the coefficients. +during each solve. The type of the result follows the evaluation point, not the basis or the coefficients. This holds +for the buffer forms and for the allocating forms alike: `evaluate`, `evaluate_all` and `Spline` +evaluation on a one-dimensional basis, `TensorProductBasis`, `RecombinedBSplineBasis` and +`PolarSplineBasis`, including the forms for one basis function and the forms for a vector of +points. A `Float64` basis with `Float64` coefficients at a `Float32` point returns a `Float32` from +`evaluate` and from `evaluate!`. `evaluate_all` returns buffers of `Float32`. Before this release +the allocating forms returned a `Float64` there, because they took the promotion of the basis and +the point. A `ForwardDiff.Dual` point returns a `ForwardDiff.Dual`. Complex coefficients make the +result complex in the type of the point: a `ComplexF64` coefficient array at a `Float32` point +returns a `ComplexF32`. A point that is not a floating-point type, such as an `Int` or a `Rational`, +cannot hold the result. For such a point the promotion of the basis type and the point type stands: +an `Int` or a `Rational` point on a `Float64` basis returns a `Float64`, and an `Int` point on a +`Float32` basis returns a `Float32`. This keeps the result of every `Int` point as it was, and needs +no new rule for them. A `Float64` point is unchanged: same type, same value, same bytes and speed. +Measured, the `@allocated` bytes are identical before and after for the one-dimensional, +tensor-product, recombined and polar forms, and the timings are equal within noise. `evaluate!` computes in the buffer's type and returns it. The buffer must hold `typeof(x)` (`eltype(x)` for the tensor-product form), or an `ArgumentError` is raised. So a `Float64` basis with `Float64` coefficients at a `Float32` point works in a `Float32` buffer. It returns a diff --git a/src/bspline.jl b/src/bspline.jl index c38b087..72190e7 100644 --- a/src/bspline.jl +++ b/src/bspline.jl @@ -1,6 +1,19 @@ @inline _evaltype(::Type{T}, ::Type{S}) where {T, S} = promote_type(T, S) +# The type of the value of an evaluation at a point of type `P` with coefficients, or a basis, of +# type `C`. It follows the point: a `Float32` point gives a `Float32` on a `Float64` basis, and +# a `ForwardDiff.Dual` point a `Dual`; complex coefficients make it complex in that type. A point +# that is not a floating-point type, an `Int` or a `Rational`, cannot hold the result, so the +# promotion of `C` and `P` stands. +@inline function _resulttype(::Type{C}, ::Type{P}) where {C, P} + if P <: Real && float(P) === P + C <: Complex ? Complex{P} : P + else + promote_type(C, P) + end +end + # Whether `x` sits at the value the topmost knot span is to be closed at. `nothing` is the # answer for a basis that has no such point -- a periodic one, whose knot vector extends # beyond the domain in both directions -- and it is a type, not a value, so the comparison @@ -564,6 +577,13 @@ The `d`-th derivative of the `j`-th basis function of `b` at `x`, or of the spli `x` may also be a vector, in which case a vector is returned. +The type of the result follows the point, as it does for [`evaluate!`](@ref): a `Float64` basis +at a `Float32` point returns a `Float32`, and a `ForwardDiff.Dual` point returns a `Dual`. +Complex coefficients make the result complex in that type. A point that is not a +floating-point type, such as an `Int` or a `Rational`, cannot hold the result, so the +promotion of the type of the basis and the type of the point stands: an `Int` point on a +`Float64` basis returns a `Float64`. + The explicit derivative order is what the lazy `b'` products cannot give: the discrete brackets need `d` up to `3`, and stacking `Derivative` that deep does not compose. @@ -593,7 +613,8 @@ two are within a third of each other at ``p = 3`` and a factor of 150 apart at ` function evaluate(b::BSplineBasis{T}, j::Integer, x::Number, d::Integer = 0) where {T} @boundscheck (1 ≤ j ≤ nbasis(b)) || throw(BoundsError(b, j)) d ≥ 0 || throw(ArgumentError("the derivative order must be non-negative, got d = $(d)")) - _bspline(b.knots, Int(j), b.p, x, Int(d), _closeat(b)) + convert(_resulttype(T, typeof(x)), _bspline( + b.knots, Int(j), b.p, x, Int(d), _closeat(b))) end function evaluate(b::PeriodicBSplineBasis{T}, j::Integer, x::Number, @@ -616,7 +637,7 @@ function evaluate(b::PeriodicBSplineBasis{T}, j::Integer, x::Number, for shift in (-L, zero(L), L) v += _bspline(b.knots, k, b.p, x̃ + shift, Int(d), nothing) end - return v + return convert(_resulttype(T, typeof(x)), v) end function evaluate(b::AbstractBSplineBasis, j::Integer, X::AbstractVector, d::Integer = 0) @@ -625,7 +646,7 @@ end function evaluate(b::AbstractBSplineBasis{T}, û::AbstractVector{S}, x::Number, d::Integer = 0) where {T, S} - R = _evaltype(promote_type(T, S), typeof(x)) + R = _resulttype(promote_type(T, S), typeof(x)) evaluate!(Vector{R}(undef, local_width(b)), b, û, x, d) end @@ -698,7 +719,7 @@ function evaluate(b::AbstractBSplineBasis{T}, û::AbstractVector{S}, X::Abstract d::Integer = 0) where {T, S} length(û) == nbasis(b) || throw(DimensionMismatch( "the coefficient vector has $(length(û)) entries but the basis has $(nbasis(b))")) - R = _evaltype(promote_type(T, S), eltype(X)) + R = _resulttype(promote_type(T, S), eltype(X)) # One buffer for the whole sweep. Evaluating point by point through the scalar method # takes one per point, which for a plot or a diagnostic over a fine grid is the entire @@ -925,14 +946,16 @@ end evaluate_all(b::AbstractBSplineBasis, x, d = 0) The allocating form of [`evaluate_all!`](@ref): returns `(j₀, values)` with a freshly allocated -buffer of [`local_width`](@ref) entries. +buffer of [`local_width`](@ref) entries. The element type of the buffer is the type of the +point, `typeof(x)`, unless that is not a floating-point type: then it is the promotion of the +type of the basis and the type of the point. Convenient at a call site that runs once. A particle loop should hold its own buffer and call [`evaluate_all!`](@ref), which allocates nothing, or [`evaluate!`](@ref) for the value of a spline at the point. """ function evaluate_all(b::AbstractBSplineBasis{T}, x::Number, d::Integer = 0) where {T} - R = _evaltype(T, typeof(x)) + R = _resulttype(T, typeof(x)) # `local_width`, not `degree(b) + 1`: a recombined basis has a wider block at the two # ends, because a recombined function there spans the union of two parent supports. values = Vector{R}(undef, local_width(b)) diff --git a/src/polar.jl b/src/polar.jl index d0c6b61..981f781 100644 --- a/src/polar.jl +++ b/src/polar.jl @@ -450,6 +450,12 @@ end The mixed derivative ``\partial_s^{d_1} \partial_\theta^{d_2}`` at `x = (s, θ)` of the `k`-th basis function, or of the spline with coefficient vector `û`. +The type of the result follows the point, `eltype(x)`, as it does for [`evaluate!`](@ref): a +`Float64` basis at a `Float32` point returns a `Float32`. Complex coefficients make the result +complex in that type. A point that is not a floating-point type, such as an `Int` or a +`Rational`, cannot hold the result, so the promotion of the type of the basis and the type of +the point stands. + ```jldoctest julia> B = PolarSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3), PeriodicBSplineBasis(UniformMesh(16, 0 .. 2π), 3)); @@ -479,7 +485,7 @@ function evaluate(B::PolarSplineBasis{T}, k::Integer, x, for t in nzrange(R, Int(k)) v += vals[t] * evaluate(B.parent, rows[t], x, d) end - return v + return convert(_resulttype(T, eltype(x)), v) end function evaluate(B::PolarSplineBasis{T}, û::AbstractVector{S}, x, @@ -504,7 +510,9 @@ Base.getindex(B::PolarSplineBasis, x, k) = evaluate(B, k, x) evaluate_all(B::PolarSplineBasis, x, d = (0, 0)) The basis functions that are nonzero at `x = (s, θ)`, as a pair `(idx, values)` of the global -indices and the matching values of ``\partial_s^{d_1} \partial_\theta^{d_2} \Psi``. +indices and the matching values of ``\partial_s^{d_1} \partial_\theta^{d_2} \Psi``. The element +type of `values` is the type of the point, `eltype(x)`, unless that is not a floating-point +type: then it is the promotion of the type of the basis and the type of the point. Unlike the tensor-product and one-dimensional forms, which return the **first index** of a contiguous block, this returns the indices themselves. The polar index set is not a product @@ -530,7 +538,10 @@ julia> length(idx), sum(vals) ≈ 1 function evaluate_all(B::PolarSplineBasis{T}, x, d::NTuple{2, Int} = (0, 0)) where {T} radial, angular = bases(B.parent) Ns = nbasis(radial) - R = _evaltype(T, eltype(x)) + R = _resulttype(T, eltype(x)) + # The pole sums mix the blocks with the `T` entries of `λ`, so they run in the promoted type + # and `push!` converts each value to `R`. + V = _evaltype(T, eltype(x)) i₀, bs = evaluate_all(radial, x[1], d[1]) j₀, bθ = evaluate_all(angular, x[2], d[2]) @@ -548,14 +559,14 @@ function evaluate_all(B::PolarSplineBasis{T}, x, d::NTuple{2, Int} = (0, 0)) whe # rather than over all of N_θ, which is what keeps this O(p²) rather than O(N_θ). if i₀ ≤ 2 for k in 1:3 - v = zero(R) + v = zero(V) for ti in eachindex(bs) i = basis_index(radial, i₀ + ti - 1) (i == 1 || i == 2) || continue - aθ = zero(R) + aθ = zero(V) for tj in eachindex(bθ) j = basis_index(angular, j₀ + tj - 1) - aθ += bθ[tj] * (i == 1 ? one(R) / 3 : B.λ[k, j]) + aθ += bθ[tj] * (i == 1 ? one(V) / 3 : B.λ[k, j]) end v += bs[ti] * aθ end diff --git a/src/recombined.jl b/src/recombined.jl index 4801da9..3d67443 100644 --- a/src/recombined.jl +++ b/src/recombined.jl @@ -290,7 +290,7 @@ function evaluate(b::RecombinedBSplineBasis{T}, j::Integer, x::Number, for t in nzrange(R, Int(j)) v += vals[t] * evaluate(b.parent, rows[t], x, d) end - return v + return convert(_resulttype(T, typeof(x)), v) end # The Greville abscissa of a recombined function is taken to be that of the parent function it diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 3da328d..0752526 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -178,6 +178,12 @@ indexed by `I`, or of the spline with coefficient array `û`. indexable `D`-vector. `d` is a per-axis tuple of derivative orders; the gradient component ``\partial_k`` is `d = ntuple(i -> i == k ? 1 : 0, D)`. +The type of the result follows the point, `eltype(x)`, as it does for [`evaluate!`](@ref): a +`Float64` basis at a `Float32` point returns a `Float32`, and a `ForwardDiff.Dual` point returns +a `Dual`. Complex coefficients make the result complex in that type. A point that is not a +floating-point type, such as an `Int` or a `Rational`, cannot hold the result, so the +promotion of the type of the basis and the type of the point stands. + ```jldoctest julia> B = BSplineBasis(UniformMesh(8, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(8, 0 .. 1), 2); @@ -212,7 +218,7 @@ end function evaluate(B::TensorProductBasis{T, D}, û::AbstractArray{S, D}, x, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {T, D, S} - R = _evaltype(promote_type(T, S), eltype(x)) + R = _resulttype(promote_type(T, S), eltype(x)) bufs = ntuple(k -> Vector{R}(undef, local_width(B.bases[k])), D) evaluate!(bufs, B, û, x, d) end @@ -312,7 +318,7 @@ end function evaluate_all(B::TensorProductBasis{T, D}, x, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {T, D} - R = _evaltype(T, eltype(x)) + R = _resulttype(T, eltype(x)) bufs = ntuple(k -> Vector{R}(undef, local_width(B.bases[k])), D) j₀ = evaluate_all!(bufs, B, x, d) return (j₀, bufs) diff --git a/test/integration/workspaces.jl b/test/integration/workspaces.jl index 8a2f2c0..96410ec 100644 --- a/test/integration/workspaces.jl +++ b/test/integration/workspaces.jl @@ -188,6 +188,105 @@ end end end +@testset "$(rpad("the allocating forms return the type of the point",76))" begin + # The reference is the Float64 evaluation at the same point, rounded to Float32. The sum + # cancels, so the absolute floor is the size of its terms: the coefficients, times the size + # of a basis function after d derivatives on 8 cells. + for b in bases_1d(Float64) + û = randn(nbasis(b)) + s = Spline(b, û) + j = nbasis(b) ÷ 2 + X = Float32[0.2, 0.37, 0.9] + for x in Float32[0.37, 0.5, 0, 1, -0.25, 1.75], d in 0:2 + + tol = 8 * eps(Float32) * (8 * degree(b))^d + near(v, ref) = isapprox(v, ref; rtol = 8 * eps(Float32), + atol = tol * max(1, maximum(abs, û))) + v = @inferred evaluate(b, û, x, d) + @test v isa Float32 + @test near(v, Float32(evaluate(b, û, Float64(x), d))) + @test near(v, evaluate!(zeros(Float32, local_width(b)), b, û, x, d)) + w = @inferred evaluate(s, x, d) + @test w isa Float32 + @test w == v + vj = @inferred evaluate(b, j, x, d) + @test vj isa Float32 + @test near(vj, Float32(evaluate(b, j, Float64(x), d))) + j₀, vals = @inferred evaluate_all(b, x, d) + @test vals isa Vector{Float32} + @test all(near(vals[t], Float32(evaluate_all(b, Float64(x), d)[2][t])) + for t in eachindex(vals)) + end + for d in 0:1 + @test (@inferred evaluate(b, û, X, d)) isa Vector{Float32} + @test (@inferred evaluate(b, j, X, d)) isa Vector{Float32} + @test evaluate(b, û, X, d) == [evaluate(b, û, x, d) for x in X] + end + + # a Float64 point keeps Float64, and so does a Float64 basis at a Float64 point + @test (@inferred evaluate(b, û, 0.37)) isa Float64 + @test (@inferred evaluate(b, j, 0.37)) isa Float64 + @test (@inferred evaluate(b, û, [0.2, 0.37])) isa Vector{Float64} + @test evaluate_all(b, 0.37)[2] isa Vector{Float64} + + # Complex coefficients: the result is complex in the type of the point + ûc = û .+ im .* randn(nbasis(b)) + vc = @inferred evaluate(b, ûc, 0.37f0) + @test vc isa ComplexF32 + @test isapprox(vc, ComplexF32(evaluate(b, ûc, 0.37)); rtol = 8 * eps(Float32), + atol = 8 * eps(Float32) * maximum(abs, ûc)) + @test (@inferred evaluate(b, ûc, 0.37)) isa ComplexF64 + @test (@inferred evaluate(b, ûc, dual(0.37))) isa Complex{Dual12} + + # A Dual point keeps its type + @test (@inferred evaluate(b, û, dual(0.37))) isa Dual12 + @test (@inferred evaluate(b, j, dual(0.37))) isa Dual12 + @test (@inferred evaluate_all(b, dual(0.37)))[2] isa Vector{Dual12} + + # A point that is not a floating-point type cannot hold the result: the promotion of the + # basis and the point stands, so a Float64 basis returns Float64 + @test (@inferred evaluate(b, û, 1)) isa Float64 + @test (@inferred evaluate(b, û, 1 // 3)) isa Float64 + @test (@inferred evaluate(b, j, 1)) isa Float64 + end + + # a Float32 basis keeps its type at an Int point, and a Float64 point stays Float64 + b32 = bases_1d(Float32)[2] + û32 = randn(Float32, nbasis(b32)) + @test (@inferred evaluate(b32, û32, 1)) isa Float32 + @test (@inferred evaluate(b32, û32, 0.37)) isa Float64 + @test (@inferred evaluate(b32, û32, 0.37f0)) isa Float32 + + for B in (basis_2d(Float64), basis_3d(Float64)) + D = ndims(B) + û = randn(size(B)...) + x = D == 2 ? (0.37f0, 0.4f0) : (0.37f0, 1.2f0, -0.3f0) + x64 = Float64.(x) + for d in (ntuple(_ -> 0, D), ntuple(k -> k == 1 ? 1 : 0, D)) + atol = 8 * eps(Float32) * maximum(abs, û) * 24^sum(d) + v = @inferred evaluate(B, û, x, d) + @test v isa Float32 + @test isapprox(v, Float32(evaluate(B, û, x64, d)); rtol = 8 * eps(Float32), atol) + w = @inferred evaluate(Spline(B, û), x, d) + @test w isa Float32 + @test w == v + I = ntuple(_ -> 2, D) + vi = @inferred evaluate(B, I, x, d) + @test vi isa Float32 + @test isapprox(vi, Float32(evaluate(B, I, x64, d)); rtol = 8 * eps(Float32), + atol = 8 * eps(Float32) * 24^sum(d)) + j₀, bufs = @inferred evaluate_all(B, x, d) + @test all(buf -> buf isa Vector{Float32}, bufs) + end + @test (@inferred evaluate(B, û, x64)) isa Float64 + @test (@inferred evaluate_all(B, x64))[2][1] isa Vector{Float64} + ûc = û .+ im .* randn(size(B)...) + @test (@inferred evaluate(B, ûc, x)) isa ComplexF32 + @test (@inferred evaluate(B, û, map(dual, x))) isa Dual12 + @test (@inferred evaluate(B, û, map(Int ∘ round, x))) isa Float64 + end +end + @testset "$(rpad("evaluate! on hostile arguments matches evaluate",76))" begin for T in (Float32, Float64), b in bases_1d(T) diff --git a/test/polar.jl b/test/polar.jl index 436490f..a5d6994 100644 --- a/test/polar.jl +++ b/test/polar.jl @@ -403,3 +403,49 @@ Random.seed!(0x2f7a91c4) 0.5 end end + +@testset "$(rpad("the allocating forms return the type of the point",80))" begin + radial = BSplineBasis(UniformMesh(10, 0 .. 1), 3) + angular = PeriodicBSplineBasis(UniformMesh(16, 0 .. 2π), 3) + B = PolarSplineBasis(radial, angular) + û = randn(nbasis(B)) + atol(d) = 8 * eps(Float32) * 24^sum(d) + + for x in ((0.0f0, 0.3f0), (0.01f0, 1.1f0), (0.5f0, 2.3f0)), + d in ((0, 0), (1, 0), (0, 1)) + + x64 = Float64.(x) + for k in (1, 3, 5, nbasis(B)) + v = @inferred evaluate(B, k, x, d) + @test v isa Float32 + @test isapprox(v, Float32(evaluate(B, k, x64, d)); rtol = 8 * eps(Float32), + atol = atol(d)) + end + v = @inferred evaluate(B, û, x, d) + @test v isa Float32 + @test isapprox(v, Float32(evaluate(B, û, x64, d)); rtol = 8 * eps(Float32), + atol = maximum(abs, û) * atol(d)) + idx, vals = @inferred evaluate_all(B, x, d) + idx64, vals64 = evaluate_all(B, x64, d) + @test vals isa Vector{Float32} + @test idx == idx64 + @test all(isapprox(vals[t], Float32(vals64[t]); rtol = 8 * eps(Float32), + atol = atol(d)) for t in eachindex(vals)) + end + @test (@inferred evaluate(B, û, [(0.5f0, 2.3f0), (0.2f0, 1.0f0)])) isa Vector{Float32} + + # a Float64 point keeps Float64 + @test (@inferred evaluate(B, 5, (0.5, 2.3))) isa Float64 + @test (@inferred evaluate(B, û, (0.5, 2.3))) isa Float64 + @test (@inferred evaluate_all(B, (0.5, 2.3)))[2] isa Vector{Float64} + + # complex coefficients, in the type of the point + ûc = û .+ im .* randn(nbasis(B)) + vc = @inferred evaluate(B, ûc, (0.5f0, 2.3f0)) + @test vc isa ComplexF32 + @test isapprox(vc, ComplexF32(evaluate(B, ûc, (0.5, 2.3))); rtol = 8 * eps(Float32), + atol = 8 * eps(Float32) * maximum(abs, ûc)) + + # a point that is not a floating-point type: the promotion of the basis and the point + @test (@inferred evaluate(B, û, (1 // 2, 1 // 3))) isa Float64 +end From 3b2f3f7b635e5e919465de3f778f33546fdb2497 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 13:50:37 +0100 Subject: [PATCH 10/12] Review threads: exact knot differences in narrow buffers, floating-point buffers only, K7 and K8 scope Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 37 ++++++++++++++++++----------- src/bspline.jl | 28 +++++++++++++--------- src/tensorproduct.jl | 8 +++++-- test/integration/workspaces.jl | 43 ++++++++++++++++++++++++++++++++++ 4 files changed, 89 insertions(+), 27 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d8ee1f5..fde5e2b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,12 +40,18 @@ no new rule for them. A `Float64` point is unchanged: same type, same value, sam Measured, the `@allocated` bytes are identical before and after for the one-dimensional, tensor-product, recombined and polar forms, and the timings are equal within noise. `evaluate!` computes in the buffer's type and returns it. The buffer must hold `typeof(x)` -(`eltype(x)` for the tensor-product form), or an `ArgumentError` is raised. So a `Float64` basis -with `Float64` coefficients at a `Float32` point works in a `Float32` buffer. It returns a -`Float32` and allocates 0 bytes, measured with `@allocated` after a warm call on a cubic clamped -basis. A `Float64` point into a `Float32` buffer raises, and so does a `ForwardDiff.Dual` point -into a `Float64` buffer. The allocating forms stay, -and build their buffers and call these. +(`eltype(x)` for the tensor-product form), or an `ArgumentError` is raised. It must also have a +floating-point element type, because the knot arithmetic is fractional: an `Int` buffer raises an +`ArgumentError` even at an `Int` point, where the allocating form returns a `Float64`. So a +`Float64` basis with `Float64` coefficients at a `Float32` point works in a `Float32` buffer. It +returns a `Float32` and allocates 0 bytes, measured with `@allocated` after a warm call on a cubic +clamped basis. A `Float64` point into a `Float32` buffer raises, and so does a `ForwardDiff.Dual` +point into a `Float64` buffer. The allocating forms stay, and build their buffers and call these. +`evaluate_all!` forms each difference of a knot and the point in the promotion of their own types, +and converts only the difference to the buffer type. A `Float64` basis on +`UniformMesh(8, 1e8 .. (1e8 + 8))` at the `Float32` point `1e8` into a `Float32` buffer therefore +returns the value of the span. The knots `1e8` and `1e8 + 1` are the same `Float32`, so converting +them before the subtraction divides by zero and returns `NaN`. Bytes per call, measured with `@allocated` after a warm call through a function barrier, in a cold process on Julia 1.13.3, on `cb82e6f` and on this release. The one-dimensional basis is @@ -68,21 +74,24 @@ before, because each builds two buffers of the largest intermediate array, and ` also the buffers of the mass solve. `mass_solve!` on a `KroneckerMass` checks the lengths of `y` and `x` before it copies `x`, and raises a `DimensionMismatch` for a `y` or an `x` whose length is not that of the operator. Where the element types of the arguments are equal, the results of -the allocating forms are bit-for-bit those of `cb82e6f`. The three-argument `l2_projection!` now +the allocating forms are bit-for-bit those of `cb82e6f`. The three-argument `l2_projection!` contracts in `eltype(û)`, not in the promoted type of the quadrature and the sample, so where -these types differ its result can differ from that of `cb82e6f` by the round-off of the narrower -type. +these types differ and the call runs, its result can differ from that of `cb82e6f` by the round-off of the narrower type. A `Float32` `û` does not run on a +`Float64` quadrature with a periodic axis on a `UniformMesh`, in `cb82e6f` and in this release: the +`CirculantMass` factor raises a `MethodError` in its `Float64` FFT plan (known issue K7). A periodic +axis on a `GradedMesh` runs. `ForwardDiff` joins the test dependencies: the tests evaluate with `ForwardDiff.Dual` arguments and buffers, which also allocate nothing. ### Bug Fixes -`mass_solve!` on a `KroneckerMass`, in the three-argument form and so in `ldiv!` and `\`, accepted -an `x` shorter than `y` in 0.3.4, with no error. It copied `x` into the start of `y` and solved the -whole of `y`, so the solve ran on whatever the rest of `y` held. It now raises a -`DimensionMismatch` for an `x` whose length is not that of `y`. This fixes known issue K8. Earlier -releases were not checked. +`mass_solve!(y, op, x)` on a `KroneckerMass`, and so `ldiv!(y, op, x)`, accepted an `x` shorter +than `y` in 0.3.4, with no error. It copied `x` into the start of `y` and solved the whole of `y`, +so the solve ran on whatever the rest of `y` held. It raises a `DimensionMismatch` for an `x` whose +length is not that of `y`. `\` is not affected: it builds its output with `similar(x)`, and a +wrong-length `x` raises a `DimensionMismatch` there in 0.3.4 and in this release. This fixes known +issue K8. Earlier releases were not checked. `mass_operator` of a periodic basis on a `UniformMesh`, and `CirculantMass` itself, raised `MethodError: … plan_rfft … is ambiguous` whenever FFTA 0.3.2 or later was loaded beside FFTW. diff --git a/src/bspline.jl b/src/bspline.jl index 72190e7..1cbf7bb 100644 --- a/src/bspline.jl +++ b/src/bspline.jl @@ -661,9 +661,11 @@ caller's buffer `buf`, which allocates nothing. contents on entry do not matter, and on return they are the local block of [`evaluate_all!`](@ref). Its element type is the type the sum is computed in, and the type of the result. It follows the point, not the basis or the coefficients: it must hold `typeof(x)`, -or an `ArgumentError` is raised. A `Float64` basis with `Float64` coefficients at a `Float32` -point takes a `Float32` buffer and returns a `Float32`. A `ForwardDiff.Dual` point takes a -buffer of `Dual`s, and a `Float64` point does not take a `Float32` buffer. +or an `ArgumentError` is raised. It must also have a floating-point element type, because the +knot arithmetic is fractional: an `Int` buffer raises an `ArgumentError` even at an `Int` point. +A `Float64` basis with `Float64` coefficients at a `Float32` point takes a `Float32` buffer and +returns a `Float32`. A `ForwardDiff.Dual` point takes a buffer of `Dual`s, and a `Float64` point +does not take a `Float32` buffer. A buffer belongs to one task. Two tasks that evaluate at once need a buffer each. @@ -686,6 +688,9 @@ function evaluate!(buf::AbstractVector{B}, b::AbstractBSplineBasis{T}, "see `local_width`")) promote_type(B, typeof(x)) === B || throw(ArgumentError( "evaluate! needs a buffer that holds the point type $(typeof(x)), got a buffer of $(B)")) + float(B) === B || throw(ArgumentError( + "evaluate! needs a buffer of a floating-point type, because the knot arithmetic is " * + "fractional, got a buffer of $(B)")) _evaluate_block!(buf, b, û, x, d) end @@ -892,13 +897,14 @@ function evaluate_all!(values::AbstractVector{R}, b::AbstractBSplineBasis, x::Nu q = p - Int(d) - # The arithmetic runs in the real type of the element type of `values`: the cell is found - # above in the type of the basis, and the knots and the point are converted once the cell - # is known. A `Float64` basis into a `Float32` buffer therefore does a `Float32` sum, with - # no `Float64` intermediate to convert back. `real`, so that a complex buffer keeps real - # knots and the comparisons of the derivative step stay real. + # The recursion runs in the real type of the element type of `values`, but every knot + # difference is formed first, in the promotion of the type of the knots and the type of the + # point, and only the difference is converted. A `Float64` basis into a `Float32` buffer + # therefore does a `Float32` sum on exact distances: converting the knots and the point + # first would round two knots far from the origin to one `Float32` and divide by zero. + # `real`, so that a complex buffer keeps real knots and the comparisons of the derivative + # step stay real. Q = real(R) - x̃ = convert(Q, x̃) # De Boor's triangular scheme for the q+1 nonzero splines of degree q on span k. # values[1:q+1] holds φ_{k-q}^q .. φ_k^q afterwards. @@ -908,8 +914,8 @@ function evaluate_all!(values::AbstractVector{R}, b::AbstractBSplineBasis, x::Nu for r in 1:j # right = kv[k+r] - x, left = kv[k+1-(j-r+1)] - ... written out below to keep # the index arithmetic of the standard form visible. - right = convert(Q, kv[k + r]) - x̃ - left = x̃ - convert(Q, kv[k + 1 - j + r - 1]) + right = convert(Q, kv[k + r] - x̃) + left = convert(Q, x̃ - kv[k + 1 - j + r - 1]) temp = values[r] / (right + left) values[r] = saved + right * temp saved = left * temp diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 0752526..099b8b3 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -234,8 +234,9 @@ which allocates nothing. their contents on entry do not matter, and on return they are the blocks of [`evaluate_all!`](@ref). Their element type is the type the sum is computed in, and the type of the result. It follows the point, not the basis or the coefficients: it must hold `eltype(x)`, -or an `ArgumentError` is raised. A `Float64` basis at a `Float32` point takes `Float32` -buffers and returns a `Float32`. +or an `ArgumentError` is raised. It must also have a floating-point element type, because the +knot arithmetic is fractional: `Int` buffers raise an `ArgumentError` even at an `Int` point. A +`Float64` basis at a `Float32` point takes `Float32` buffers and returns a `Float32`. The buffers belong to one task. Two tasks that evaluate at once need a tuple each. """ @@ -249,6 +250,9 @@ function evaluate!(bufs::NTuple{D, <:AbstractVector}, B::TensorProductBasis{T, D R = eltype(first(bufs)) promote_type(R, eltype(x)) === R || throw(ArgumentError( "evaluate! needs buffers that hold the point type $(eltype(x)), got buffers of $(R)")) + float(R) === R || throw(ArgumentError( + "evaluate! needs buffers of a floating-point type, because the knot arithmetic is " * + "fractional, got buffers of $(R)")) # Only the local block contributes, so this is the fast path rather than a sum over the # whole array. `evaluate_all!` gives the p_d+1 nonzero factors on each axis; their outer diff --git a/test/integration/workspaces.jl b/test/integration/workspaces.jl index 96410ec..fca135d 100644 --- a/test/integration/workspaces.jl +++ b/test/integration/workspaces.jl @@ -137,6 +137,41 @@ end @test startswith(message(() -> evaluate!(zeros(w), b, û, dual(0.37))), "ArgumentError: evaluate! needs a buffer that holds the point type $(Dual12), " * "got a buffer of Float64") + # an integer buffer cannot hold the fractional knot arithmetic, even at an Int point, + # where the allocating form returns a Float64 + @test_throws ArgumentError evaluate!(zeros(Int, w), b, û, 0) + @test startswith(message(() -> evaluate!(zeros(Int, w), b, û, 0)), + "ArgumentError: evaluate! needs a buffer of a floating-point type, because the " * + "knot arithmetic is fractional, got a buffer of Int64") + @test_throws ArgumentError evaluate!(zeros(Int, w), Spline(b, û), 0) + @test evaluate(b, û, 0) isa Float64 + # an Int point does take a floating-point buffer + @test evaluate!(zeros(w), b, û, 0) == evaluate(b, û, 0) + end +end + +@testset "$(rpad("evaluate! into a narrow buffer subtracts the knots before narrowing",76))" begin + # Two knots 1 apart at 1e8 are one Float32. The distances are exact in the Float64 basis, + # so the Float32 result is the Float64 one, rounded. + b = BSplineBasis(UniformMesh(8, 1.0e8 .. (1.0e8 + 8)), 3) + û = collect(1.0:nbasis(b)) + buf = zeros(Float32, local_width(b)) + for x in (Float32(1.0e8), Float32(1.0e8 + 8)), d in 0:2 + + v = @inferred evaluate!(buf, b, û, x, d) + @test v isa Float32 + @test isfinite(v) + @test isapprox(v, Float32(evaluate(b, û, Float64(x), d)); rtol = 8 * eps(Float32), + atol = 8 * eps(Float32) * maximum(abs, û)) + @test isfinite(evaluate(b, û, x, d)) + end + # an interior point, where the Float32 grid is fine enough to hold it + b = BSplineBasis(UniformMesh(8, 1.0e6 .. (1.0e6 + 8)), 3) + x = Float32(1.0e6 + 3.5) + for d in 0:2 + v = evaluate!(zeros(Float32, local_width(b)), b, û, x, d) + @test isapprox(v, Float32(evaluate(b, û, Float64(x), d)); rtol = 8 * eps(Float32), + atol = 8 * eps(Float32) * maximum(abs, û) * 8^d) end end @@ -373,6 +408,14 @@ end @test startswith(message(() -> evaluate!(narrow, B, û, first(points))), "ArgumentError: evaluate! needs buffers that hold the point type Float64, " * "got buffers of Float32") + # integer buffers, which cannot hold the fractional knot arithmetic, even at an + # Int point + ints = map(w -> zeros(Int, w), local_width(B)) + xi = ntuple(_ -> 0, D) + @test_throws ArgumentError evaluate!(ints, B, û, xi) + @test startswith(message(() -> evaluate!(ints, B, û, xi)), + "ArgumentError: evaluate! needs buffers of a floating-point type, because " * + "the knot arithmetic is fractional, got buffers of Int64") end end end From 037f9c3e7bf4616c552e56f03b1de0ee3adddd7e Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 17:45:32 +0100 Subject: [PATCH 11/12] Narrow the result type only for IEEE floats; 0.4.0 breaking changes The result of an evaluation takes the precision of the point only where the basis and the coefficients are IEEE floats or complex numbers of them and the point is an IEEE float. Every other type is kept by promotion, so Dual coefficients at a Float64 point return a Dual again (ForwardDiff.gradient with respect to the coefficients works), and BigFloat coefficients or a BigFloat basis return a BigFloat. The rule is written as dispatch, which the compiler resolves; a subtype test in the body allocated on every call. A recombined basis sums its parent's values in the promoted type, so a buffer wider than the point keeps the precision of the basis. Tests: Dual-coefficient gradients in 1D and in the tensor product, BigFloat and complex BigFloat coefficients, a BigFloat basis, and complex coefficients at d = 1, 2 (mutant M11). CHANGELOG: the section targets 0.4.0, with a Breaking Changes section for the narrowed result type, the DimensionMismatch of mass_solve! and ldiv!, and the buffers the new forms refuse; the K8 reference goes. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 78 +++++++++++++++++++-------------- src/bspline.jl | 68 +++++++++++++++++------------ src/polar.jl | 14 +++--- src/recombined.jl | 2 +- src/tensorproduct.jl | 20 +++++---- test/integration/workspaces.jl | 79 +++++++++++++++++++++++++++++++++- 6 files changed, 183 insertions(+), 78 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a02500..0dbf1f1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,42 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] — targeting 0.4.0 +### Breaking Changes + +**An evaluation at a narrower IEEE float point returns the precision of the point.** Where the +basis and the coefficients are IEEE floats (`Float16`, `Float32`, `Float64`) or complex numbers of +them, and the point is an IEEE float, the result takes the precision of the point. A `Float64` +basis with `Float64` coefficients at a `Float32` point returns a `Float32`, where 0.3.5 returned a +`Float64`; complex `Float64` coefficients there return a `ComplexF32`, where 0.3.5 returned a +`ComplexF64`; and `evaluate_all` returns buffers of `Float32`. This holds for `evaluate`, +`evaluate_all` and `Spline` evaluation on a one-dimensional basis, `TensorProductBasis`, +`RecombinedBSplineBasis` and `PolarSplineBasis`, including the forms for one basis function and +the forms for a vector of points. A caller that needs the 0.3.5 result passes the point in the +precision of the basis: `Float64(x)`, or `Float64.(x)` for a tuple or a vector of coordinates. +Every other type is kept by promotion of the types of the basis, the coefficients and the point, +as in 0.3.5: a `ForwardDiff.Dual` point or `Dual` coefficients return a `Dual`, so +`ForwardDiff.gradient` with respect to the coefficients works; a `BigFloat` basis or `BigFloat` +coefficients return a `BigFloat`; and an `Int` or a `Rational` point on a `Float64` basis returns +a `Float64`. A point of the precision of the basis, or a wider one, returns the type of 0.3.5. + +**`mass_solve!(y, op, x)` on a `KroneckerMass`, and so `ldiv!(y, op, x)`, raise a +`DimensionMismatch` for an `x` whose length is not that of `y`.** In 0.3.4 and 0.3.5 an `x` +shorter than `y` raised no error: the solve copied `x` into the start of `y` and solved the whole +of `y`, so it ran on whatever the rest of `y` held. An `x` longer than `y` raised a `BoundsError`. +A caller passes an `x` of `length(y)`; a caller that catches the `BoundsError` catches the +`DimensionMismatch` instead. `\` is not affected: it builds its output with `similar(x)`, and a +wrong-length `x` raises a `DimensionMismatch` there in 0.3.5 and in this release. + +**The new caller-buffer forms refuse some buffers.** These forms are new in this release, so no +call of 0.3.5 changes. `evaluate!` raises an `ArgumentError` for a buffer that does not hold the +type of the point, or whose element type is not a floating-point type; an `Int` buffer is refused +even at an `Int` point, where the allocating form returns a `Float64`. A caller allocates the +buffer in the type of the point, `zeros(typeof(x), local_width(b))`, or in a wider floating-point +type. `mass_solve!(y, op, x, work)` and `contract!(L, q, F, work, d)` raise an `ArgumentError` for +a workspace whose element type is not that of `y` or `L`, or whose sizes do not fit the operator +or the quadrature: a caller builds the workspace with `mass_workspace(op, eltype(y))` or +`contract_workspace(q, eltype(L))`. + ### New Features The hot evaluation, solve and projection paths take buffers that the caller holds, and then @@ -23,26 +59,10 @@ does it on a `TensorProductBasis` with one buffer per axis. `mass_workspace(op, `contract!(L, q, F, work, d)` and `l2_projection!(û, q, F, work)` contract and project on a `TensorProductQuadrature`. Each buffer and workspace belongs to one task. Tasks must not share a mass operator or `TensorProductQuadrature`, because a `CirculantMass` factor writes scratch space -during each solve. The type of the result follows the evaluation point, not the basis or the coefficients. This holds -for the buffer forms and for the allocating forms alike: `evaluate`, `evaluate_all` and `Spline` -evaluation on a one-dimensional basis, `TensorProductBasis`, `RecombinedBSplineBasis` and -`PolarSplineBasis`, including the forms for one basis function and the forms for a vector of -points. A `Float64` basis with `Float64` coefficients at a `Float32` point returns a `Float32` from -`evaluate` and from `evaluate!`. `evaluate_all` returns buffers of `Float32`. Before this release -the allocating forms returned a `Float64` there, because they took the promotion of the basis and -the point. A `ForwardDiff.Dual` point returns a `ForwardDiff.Dual`. Complex coefficients make the -result complex in the type of the point: a `ComplexF64` coefficient array at a `Float32` point -returns a `ComplexF32`. A point that is not a floating-point type, such as an `Int` or a `Rational`, -cannot hold the result. For such a point the promotion of the basis type and the point type stands: -an `Int` or a `Rational` point on a `Float64` basis returns a `Float64`, and an `Int` point on a -`Float32` basis returns a `Float32`. This keeps the result of every `Int` point as it was, and needs -no new rule for them. A `Float64` point is unchanged: same type, same value, same bytes and speed. -Measured, the `@allocated` bytes are identical before and after for the one-dimensional, -tensor-product, recombined and polar forms, and the timings are equal within noise. -`evaluate!` computes in the buffer's type and returns it. The buffer must hold `typeof(x)` -(`eltype(x)` for the tensor-product form), or an `ArgumentError` is raised. It must also have a -floating-point element type, because the knot arithmetic is fractional: an `Int` buffer raises an -`ArgumentError` even at an `Int` point, where the allocating form returns a `Float64`. So a +during each solve. The type of the result follows the rule of the first breaking change above, +for the buffer forms and for the allocating forms alike. `evaluate!` computes in the buffer's type +and returns it. The buffer must hold `typeof(x)` (`eltype(x)` for the tensor-product form), and it +must have a floating-point element type; the last breaking change above gives the refusals. So a `Float64` basis with `Float64` coefficients at a `Float32` point works in a `Float32` buffer. It returns a `Float32` and allocates 0 bytes, measured with `@allocated` after a warm call on a cubic clamped basis. A `Float64` point into a `Float32` buffer raises, and so does a `ForwardDiff.Dual` @@ -76,23 +96,15 @@ and `x` before it copies `x`, and raises a `DimensionMismatch` for a `y` or an ` is not that of the operator. Where the element types of the arguments are equal, the results of the allocating forms are bit-for-bit those of `cb82e6f`. The three-argument `l2_projection!` contracts in `eltype(û)`, not in the promoted type of the quadrature and the sample, so where -these types differ and the call runs, its result can differ from that of `cb82e6f` by the round-off of the narrower type. A `Float32` `û` does not run on a -`Float64` quadrature with a periodic axis on a `UniformMesh`, in `cb82e6f` and in this release: the -`CirculantMass` factor raises a `MethodError` in its `Float64` FFT plan (known issue K7). A periodic -axis on a `GradedMesh` runs. +these types differ and the call runs, its result can differ from that of `cb82e6f` by the +round-off of the narrower type. A `Float32` `û` does not run on a `Float64` quadrature with a +periodic axis on a `UniformMesh`, in `cb82e6f` and in this release: the `CirculantMass` factor +raises a `MethodError` in its `Float64` FFT plan (known issue K7). A periodic axis on a +`GradedMesh` runs. `ForwardDiff` joins the test dependencies: the tests evaluate with `ForwardDiff.Dual` arguments and buffers, which also allocate nothing. -### Bug Fixes - -`mass_solve!(y, op, x)` on a `KroneckerMass`, and so `ldiv!(y, op, x)`, accepted an `x` shorter -than `y` in 0.3.4, with no error. It copied `x` into the start of `y` and solved the whole of `y`, -so the solve ran on whatever the rest of `y` held. It raises a `DimensionMismatch` for an `x` whose -length is not that of `y`. `\` is not affected: it builds its output with `similar(x)`, and a -wrong-length `x` raises a `DimensionMismatch` there in 0.3.4 and in this release. This fixes known -issue K8. Earlier releases were not checked. - ## [0.3.5] — 2026-11-02 ### Bug Fixes diff --git a/src/bspline.jl b/src/bspline.jl index 1cbf7bb..cbb0540 100644 --- a/src/bspline.jl +++ b/src/bspline.jl @@ -1,17 +1,22 @@ @inline _evaltype(::Type{T}, ::Type{S}) where {T, S} = promote_type(T, S) -# The type of the value of an evaluation at a point of type `P` with coefficients, or a basis, of -# type `C`. It follows the point: a `Float32` point gives a `Float32` on a `Float64` basis, and -# a `ForwardDiff.Dual` point a `Dual`; complex coefficients make it complex in that type. A point -# that is not a floating-point type, an `Int` or a `Rational`, cannot hold the result, so the -# promotion of `C` and `P` stands. -@inline function _resulttype(::Type{C}, ::Type{P}) where {C, P} - if P <: Real && float(P) === P - C <: Complex ? Complex{P} : P - else - promote_type(C, P) - end +# The type of the value of an evaluation at a point of type `P`, where `C` is the type of the +# basis or the promotion of the basis and the coefficients. Where `C` is an IEEE float +# (`Float16`, `Float32`, `Float64`) or a complex number of one, and `P` is an IEEE float, the +# result takes the precision of the point: a `Float64` basis at a `Float32` point gives a +# `Float32`, and complex `Float64` coefficients there a `ComplexF32`. Every other type of the +# basis, the coefficients or the point is kept by promotion: `ForwardDiff.Dual` coefficients or a +# `Dual` point give a `Dual`, a `BigFloat` basis or `BigFloat` coefficients a `BigFloat`, and an +# `Int` or a `Rational` point the type of the basis. The rule is written as dispatch, so that it +# resolves when the caller is compiled; a subtype test against +# `Union{Base.IEEEFloat, Complex{<:Base.IEEEFloat}}` in the body is not folded, and allocates on +# every call. +@inline _resulttype(::Type{C}, ::Type{P}) where {C, P} = promote_type(C, P) +@inline _resulttype(::Type{<:Base.IEEEFloat}, ::Type{P}) where {P <: Base.IEEEFloat} = P +@inline function _resulttype(::Type{<:Complex{<:Base.IEEEFloat}}, + ::Type{P}) where {P <: Base.IEEEFloat} + Complex{P} end # Whether `x` sits at the value the topmost knot span is to be closed at. `nothing` is the @@ -577,12 +582,13 @@ The `d`-th derivative of the `j`-th basis function of `b` at `x`, or of the spli `x` may also be a vector, in which case a vector is returned. -The type of the result follows the point, as it does for [`evaluate!`](@ref): a `Float64` basis -at a `Float32` point returns a `Float32`, and a `ForwardDiff.Dual` point returns a `Dual`. -Complex coefficients make the result complex in that type. A point that is not a -floating-point type, such as an `Int` or a `Rational`, cannot hold the result, so the -promotion of the type of the basis and the type of the point stands: an `Int` point on a -`Float64` basis returns a `Float64`. +Where the basis and the coefficients are IEEE floats (`Float16`, `Float32`, `Float64`) or +complex numbers of them, and the point is an IEEE float, the result takes the precision of the +point: a `Float64` basis at a `Float32` point returns a `Float32`, and complex `Float64` +coefficients there a `ComplexF32`. Every other type is kept by promotion of the types of the +basis, the coefficients and the point: a `ForwardDiff.Dual` point or `Dual` coefficients return +a `Dual`, a `BigFloat` basis or `BigFloat` coefficients a `BigFloat`, and an `Int` point on a +`Float64` basis a `Float64`. The explicit derivative order is what the lazy `b'` products cannot give: the discrete brackets need `d` up to `3`, and stacking `Derivative` that deep does not compose. @@ -612,9 +618,15 @@ two are within a third of each other at ``p = 3`` and a factor of 150 apart at ` """ function evaluate(b::BSplineBasis{T}, j::Integer, x::Number, d::Integer = 0) where {T} @boundscheck (1 ≤ j ≤ nbasis(b)) || throw(BoundsError(b, j)) + convert(_resulttype(T, typeof(x)), _evaluate_promoted(b, j, x, d)) +end + +# The value of `evaluate(b, j, x, d)` in the promotion of the type of the basis and the type of +# the point, before the result type narrows it. A recombined basis sums these, so that a buffer +# wider than the point keeps the precision of the basis. +function _evaluate_promoted(b::BSplineBasis, j::Integer, x::Number, d::Integer) d ≥ 0 || throw(ArgumentError("the derivative order must be non-negative, got d = $(d)")) - convert(_resulttype(T, typeof(x)), _bspline( - b.knots, Int(j), b.p, x, Int(d), _closeat(b))) + _bspline(b.knots, Int(j), b.p, x, Int(d), _closeat(b)) end function evaluate(b::PeriodicBSplineBasis{T}, j::Integer, x::Number, @@ -660,12 +672,14 @@ caller's buffer `buf`, which allocates nothing. `buf` must have [`local_width`](@ref)`(b)` entries, or a `DimensionMismatch` is raised; its contents on entry do not matter, and on return they are the local block of [`evaluate_all!`](@ref). Its element type is the type the sum is computed in, and the type of -the result. It follows the point, not the basis or the coefficients: it must hold `typeof(x)`, -or an `ArgumentError` is raised. It must also have a floating-point element type, because the -knot arithmetic is fractional: an `Int` buffer raises an `ArgumentError` even at an `Int` point. -A `Float64` basis with `Float64` coefficients at a `Float32` point takes a `Float32` buffer and -returns a `Float32`. A `ForwardDiff.Dual` point takes a buffer of `Dual`s, and a `Float64` point -does not take a `Float32` buffer. +the result; [`evaluate`](@ref) allocates a buffer of the type that its result-type rule gives. +The buffer must hold `typeof(x)`, or an `ArgumentError` is raised. It must also have a +floating-point element type, because the knot arithmetic is fractional: an `Int` buffer raises +an `ArgumentError` even at an `Int` point. A `Float64` basis with `Float64` coefficients at a +`Float32` point takes a `Float32` buffer and returns a `Float32`. A `ForwardDiff.Dual` point +takes a buffer of `Dual`s, and a `Float64` point does not take a `Float32` buffer. Coefficients +that the buffer cannot hold, such as `Dual` coefficients in a `Float64` buffer, raise the +conversion error of the sum. A buffer belongs to one task. Two tasks that evaluate at once need a buffer each. @@ -953,8 +967,8 @@ end The allocating form of [`evaluate_all!`](@ref): returns `(j₀, values)` with a freshly allocated buffer of [`local_width`](@ref) entries. The element type of the buffer is the type of the -point, `typeof(x)`, unless that is not a floating-point type: then it is the promotion of the -type of the basis and the type of the point. +point, `typeof(x)`, where the basis and the point are IEEE floats; otherwise it is the promotion +of the type of the basis and the type of the point. Convenient at a call site that runs once. A particle loop should hold its own buffer and call [`evaluate_all!`](@ref), which allocates nothing, or [`evaluate!`](@ref) for the value of a diff --git a/src/polar.jl b/src/polar.jl index 981f781..84349d8 100644 --- a/src/polar.jl +++ b/src/polar.jl @@ -450,11 +450,11 @@ end The mixed derivative ``\partial_s^{d_1} \partial_\theta^{d_2}`` at `x = (s, θ)` of the `k`-th basis function, or of the spline with coefficient vector `û`. -The type of the result follows the point, `eltype(x)`, as it does for [`evaluate!`](@ref): a -`Float64` basis at a `Float32` point returns a `Float32`. Complex coefficients make the result -complex in that type. A point that is not a floating-point type, such as an `Int` or a -`Rational`, cannot hold the result, so the promotion of the type of the basis and the type of -the point stands. +The type of the result follows the rule of the one-dimensional [`evaluate`](@ref), with +`eltype(x)` as the type of the point: an IEEE float basis with IEEE float or complex coefficients +at a `Float32` point returns a `Float32` or a `ComplexF32`, and every other type of the basis, +the coefficients or the point, such as a `ForwardDiff.Dual`, a `BigFloat` or an `Int`, is kept +by promotion. ```jldoctest julia> B = PolarSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3), @@ -511,8 +511,8 @@ Base.getindex(B::PolarSplineBasis, x, k) = evaluate(B, k, x) The basis functions that are nonzero at `x = (s, θ)`, as a pair `(idx, values)` of the global indices and the matching values of ``\partial_s^{d_1} \partial_\theta^{d_2} \Psi``. The element -type of `values` is the type of the point, `eltype(x)`, unless that is not a floating-point -type: then it is the promotion of the type of the basis and the type of the point. +type of `values` is the type of the point, `eltype(x)`, where the basis and the point are IEEE +floats; otherwise it is the promotion of the type of the basis and the type of the point. Unlike the tensor-product and one-dimensional forms, which return the **first index** of a contiguous block, this returns the indices themselves. The polar index set is not a product diff --git a/src/recombined.jl b/src/recombined.jl index 3d67443..e1b1436 100644 --- a/src/recombined.jl +++ b/src/recombined.jl @@ -358,7 +358,7 @@ function evaluate_all!(values::AbstractVector{R}, b::RecombinedBSplineBasis, x:: vals = nonzeros(Rm) for j in j₀:jn, t in nzrange(Rm, j) - values[j - j₀ + 1] += vals[t] * evaluate(b.parent, rows[t], x, d) + values[j - j₀ + 1] += vals[t] * _evaluate_promoted(b.parent, rows[t], x, d) end return j₀ end diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index 099b8b3..5f857ed 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -178,11 +178,11 @@ indexed by `I`, or of the spline with coefficient array `û`. indexable `D`-vector. `d` is a per-axis tuple of derivative orders; the gradient component ``\partial_k`` is `d = ntuple(i -> i == k ? 1 : 0, D)`. -The type of the result follows the point, `eltype(x)`, as it does for [`evaluate!`](@ref): a -`Float64` basis at a `Float32` point returns a `Float32`, and a `ForwardDiff.Dual` point returns -a `Dual`. Complex coefficients make the result complex in that type. A point that is not a -floating-point type, such as an `Int` or a `Rational`, cannot hold the result, so the -promotion of the type of the basis and the type of the point stands. +The type of the result follows the rule of the one-dimensional [`evaluate`](@ref), with +`eltype(x)` as the type of the point: an IEEE float basis with IEEE float or complex coefficients +at a `Float32` point returns a `Float32` or a `ComplexF32`, and every other type of the basis, +the coefficients or the point, such as a `ForwardDiff.Dual`, a `BigFloat` or an `Int`, is kept +by promotion. ```jldoctest julia> B = BSplineBasis(UniformMesh(8, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(8, 0 .. 1), 2); @@ -233,10 +233,12 @@ which allocates nothing. `bufs[k]` must have `local_width(bases(B)[k])` entries, or a `DimensionMismatch` is raised; their contents on entry do not matter, and on return they are the blocks of [`evaluate_all!`](@ref). Their element type is the type the sum is computed in, and the type of -the result. It follows the point, not the basis or the coefficients: it must hold `eltype(x)`, -or an `ArgumentError` is raised. It must also have a floating-point element type, because the -knot arithmetic is fractional: `Int` buffers raise an `ArgumentError` even at an `Int` point. A -`Float64` basis at a `Float32` point takes `Float32` buffers and returns a `Float32`. +the result; [`evaluate`](@ref) allocates buffers of the type that its result-type rule gives. +The buffers must hold `eltype(x)`, or an `ArgumentError` is raised. They must also have a +floating-point element type, because the knot arithmetic is fractional: `Int` buffers raise an +`ArgumentError` even at an `Int` point. A `Float64` basis at a `Float32` point takes `Float32` +buffers and returns a `Float32`. Coefficients that the buffers cannot hold, such as `Dual` +coefficients in `Float64` buffers, raise the conversion error of the sum. The buffers belong to one task. Two tasks that evaluate at once need a tuple each. """ diff --git a/test/integration/workspaces.jl b/test/integration/workspaces.jl index fca135d..bcec4c9 100644 --- a/test/integration/workspaces.jl +++ b/test/integration/workspaces.jl @@ -1,4 +1,4 @@ -using ForwardDiff: Dual +using ForwardDiff: ForwardDiff, Dual using LinearAlgebra using Random using SimpleSplines @@ -322,6 +322,83 @@ end end end +@testset "$(rpad("coefficients that are not IEEE floats keep their type",76))" begin + x = 0.37 + for b in bases_1d(Float64) + û = collect(range(0.1, 2.0; length = nbasis(b))) + φ = [evaluate(b, j, x) for j in 1:nbasis(b)] + φ′ = [evaluate(b, j, x, 1) for j in 1:nbasis(b)] + + # Dual coefficients at a Float64 point: the gradient of the value with respect to the + # coefficients is the vector of the basis values at the point. + g = ForwardDiff.gradient(c -> evaluate(b, c, x), û) + @test g isa Vector{Float64} + @test isapprox(g, φ; rtol = 8 * eps(), atol = 8 * eps()) + g′ = ForwardDiff.gradient(c -> evaluate(b, c, x, 1), û) + @test isapprox(g′, φ′; rtol = 8 * eps(), atol = 8 * eps() * 24) + @test isapprox(ForwardDiff.gradient(c -> Spline(b, c)(x), û), φ; rtol = 8 * eps(), + atol = 8 * eps()) + @test isapprox(ForwardDiff.gradient(c -> evaluate(b, c, Float32(x)), û), + [evaluate(b, j, Float32(x)) for j in 1:nbasis(b)]; rtol = 8 * eps(Float32), + atol = 8 * eps(Float32)) + + ûd = map(dual, û) + vd = @inferred evaluate(b, ûd, x) + @test vd isa Dual12 + buf = zeros(Dual12, local_width(b)) + @test (@inferred evaluate!(buf, b, ûd, x)) == vd + # A Float64 buffer cannot hold Dual coefficients: the conversion of the sum refuses them. + @test outcome(() -> evaluate!(zeros(local_width(b)), b, ûd, x)) == MethodError + + # BigFloat coefficients at a Float64 or a Float32 point return a BigFloat. + for xp in (x, Float32(x)) + vb = evaluate(b, big.(û), xp) + @test vb isa BigFloat + @test isapprox(vb, evaluate(b, û, Float64(xp)); rtol = 8 * eps()) + # Complex BigFloat coefficients return a Complex{BigFloat}. + @test evaluate(b, complex.(big.(û), big.(û)), xp) isa Complex{BigFloat} + end + end + + # A BigFloat basis at a Float64 point returns a BigFloat. + for bb in bases_1d(BigFloat) + û = big.(collect(range(0.1, 2.0; length = nbasis(bb)))) + @test evaluate(bb, û, x) isa BigFloat + @test evaluate(bb, Float64.(û), x) isa BigFloat + @test evaluate(bb, 2, x) isa BigFloat + @test evaluate_all(bb, x)[2] isa Vector{BigFloat} + end + + for B in (basis_2d(Float64), basis_3d(Float64)) + D = ndims(B) + Û = ones(size(B)...) + xB = D == 2 ? (0.3, 0.4) : (0.3, 1.2, -0.3) + g = ForwardDiff.gradient(c -> evaluate(B, c, xB), Û) + @test g isa Array{Float64, D} + @test isapprox(g, [evaluate(B, I, xB) for I in CartesianIndices(B)]; + rtol = 8 * eps(), atol = 8 * eps()) + @test evaluate(B, big.(Û), xB) isa BigFloat + end +end + +@testset "$(rpad("complex coefficients at derivative orders above zero",76))" begin + # The derivative step compares knot differences, which stay real in a complex buffer. + for b in bases_1d(Float64) + ur = randn(nbasis(b)) + ui = randn(nbasis(b)) + for x in (0.37, 0.37f0), d in 1:2 + + T = typeof(x) + v = @inferred evaluate(b, ur .+ im .* ui, x, d) + @test v isa Complex{T} + tol = 8 * eps(T) * (8 * degree(b))^d * 2 + @test isapprox(real(v), evaluate(b, ur, x, d); rtol = 8 * eps(T), atol = tol) + @test isapprox(imag(v), evaluate(b, ui, x, d); rtol = 8 * eps(T), atol = tol) + @test abs(imag(v)) > 0 || degree(b) < d + end + end +end + @testset "$(rpad("evaluate! on hostile arguments matches evaluate",76))" begin for T in (Float32, Float64), b in bases_1d(T) From 44ed23000e319c317d1966aff4af7c3b637fdde5 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 2 Nov 2026 18:02:24 +0100 Subject: [PATCH 12/12] State the result-type change by the promotion of the basis and the coefficients The Breaking Changes entry said that a point of the precision of the basis keeps the 0.3.5 type. A Float32 basis with Float64 coefficients at a Float32 point now returns Float32, where 0.3.5 returned Float64. The entry now names the change wherever the point is narrower than the promotion of the basis and the coefficients, and the remedy passes the point in that precision. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0dbf1f1..c96dcfa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,21 +9,25 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### Breaking Changes -**An evaluation at a narrower IEEE float point returns the precision of the point.** Where the -basis and the coefficients are IEEE floats (`Float16`, `Float32`, `Float64`) or complex numbers of -them, and the point is an IEEE float, the result takes the precision of the point. A `Float64` +**An evaluation at an IEEE float point returns the precision of the point.** Where the basis and +the coefficients are IEEE floats (`Float16`, `Float32`, `Float64`) or complex numbers of them, and +the point is an IEEE float, the result takes the precision of the point. This changes the result +wherever the point is narrower than the promotion of the basis and the coefficients. A `Float64` basis with `Float64` coefficients at a `Float32` point returns a `Float32`, where 0.3.5 returned a -`Float64`; complex `Float64` coefficients there return a `ComplexF32`, where 0.3.5 returned a -`ComplexF64`; and `evaluate_all` returns buffers of `Float32`. This holds for `evaluate`, +`Float64`; so does a `Float32` basis with `Float64` coefficients at a `Float32` point; complex +`Float64` coefficients there return a `ComplexF32`, where 0.3.5 returned a `ComplexF64`; and +`evaluate_all` returns buffers of `Float32`. This holds for `evaluate`, `evaluate_all` and `Spline` evaluation on a one-dimensional basis, `TensorProductBasis`, `RecombinedBSplineBasis` and `PolarSplineBasis`, including the forms for one basis function and the forms for a vector of points. A caller that needs the 0.3.5 result passes the point in the -precision of the basis: `Float64(x)`, or `Float64.(x)` for a tuple or a vector of coordinates. +precision of `promote_type(eltype(basis), eltype(û))`: for `Float64` coefficients, `Float64(x)`, or +`Float64.(x)` for a tuple or a vector of coordinates. Every other type is kept by promotion of the types of the basis, the coefficients and the point, as in 0.3.5: a `ForwardDiff.Dual` point or `Dual` coefficients return a `Dual`, so `ForwardDiff.gradient` with respect to the coefficients works; a `BigFloat` basis or `BigFloat` coefficients return a `BigFloat`; and an `Int` or a `Rational` point on a `Float64` basis returns -a `Float64`. A point of the precision of the basis, or a wider one, returns the type of 0.3.5. +a `Float64`. A point of the precision of that promotion, or a wider one, returns the type of +0.3.5. **`mass_solve!(y, op, x)` on a `KroneckerMass`, and so `ldiv!(y, op, x)`, raise a `DimensionMismatch` for an `x` whose length is not that of `y`.** In 0.3.4 and 0.3.5 an `x`