diff --git a/CHANGELOG.md b/CHANGELOG.md index 71c8bef..c96dcfa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,110 @@ 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.4.0 + +### Breaking Changes + +**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`; 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 `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 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` +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 +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. 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 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` +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 +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 | 29584 / 14912 | 0 / 0 | +| `l2_projection!(û, q, F)` | 25296 / 14640 | 29344 / 14880 | 0 / 0 | + +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 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!` +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. + +`ForwardDiff` joins the test dependencies: the tests evaluate with `ForwardDiff.Dual` arguments +and buffers, which also allocate nothing. + ## [0.3.5] — 2026-11-02 ### Bug Fixes diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 424c4ac..c4b70b1 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -71,3 +71,40 @@ `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 + +### 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` + +### 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}`, 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)`, 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/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..cbb0540 100644 --- a/src/bspline.jl +++ b/src/bspline.jl @@ -1,6 +1,24 @@ @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`, 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 # 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 +582,14 @@ 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. +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. @@ -592,6 +618,13 @@ 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)")) _bspline(b.knots, Int(j), b.p, x, Int(d), _closeat(b)) end @@ -616,7 +649,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,22 +658,66 @@ end function evaluate(b::AbstractBSplineBasis{T}, û::AbstractVector{S}, x::Number, d::Integer = 0) where {T, S} + R = _resulttype(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 computed in, and the type of +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. + +```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}, 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))")) - 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`")) + 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 # 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) @@ -652,7 +729,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 @@ -661,7 +738,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 @@ -834,6 +911,15 @@ function evaluate_all!(values::AbstractVector{R}, b::AbstractBSplineBasis, x::Nu q = p - Int(d) + # 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) + # 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) @@ -842,8 +928,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 = convert(Q, x̃ - kv[k + 1 - j + r - 1]) temp = values[r] / (right + left) values[r] = saved + right * temp saved = left * temp @@ -863,8 +949,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) @@ -880,13 +966,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)`, 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. +[`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..84349d8 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 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), 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)`, 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 @@ -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..e1b1436 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 @@ -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/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..5f857ed 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 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); @@ -212,14 +218,47 @@ end function evaluate(B::TensorProductBasis{T, D}, û::AbstractArray{S, D}, x, d::NTuple{D, Int} = ntuple(_ -> 0, D)) where {T, D, S} + 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 + +@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 computed in, and the type of +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. +""" +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))")) - 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))")) + 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 # 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 @@ -238,7 +277,7 @@ function evaluate(B::TensorProductBasis{T, D}, û::AbstractArray{S, D}, x, # 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 @@ -285,7 +324,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) @@ -380,15 +419,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 @@ -399,52 +437,120 @@ 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 +# 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 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}} + + # 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 + +@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 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 +`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. + +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)) 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)")) + length(x) == length(y) || throw(DimensionMismatch( + "the argument has $(length(x)) 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 +760,156 @@ 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)) + 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 <: Union{KroneckerMassWorkspace{T}, Nothing}} + a::Vector{T} + b::Vector{T} + mass::MW +end - # 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) +Base.eltype(::ContractWorkspace{T}) where {T} = T - # 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. +# 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 - w = quadrature_weights(q.quadratures[k]) - _scale_along!(B, w, k) + m = m ÷ Qs[k] * Ns[k] + n = max(n, m) end + return n +end + +""" + contract_workspace(q::TensorProductQuadrature, T) + +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. 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 +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 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. + +```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)) 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. 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 - Φ = 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 + _scale_along!(work.a, size(F), w[k], k) end - return reshape(B, size(q.basis)) + C = _contract_axes!(work.a, work.b, length(F), q.quadratures, d) + copyto!(L, firstindex(L), C, firstindex(C), length(L)) + return L end -function _scale_along!(A::AbstractArray, w::AbstractVector, k::Integer) - dims = size(A) +# 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 +# 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 +917,16 @@ 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, 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); @@ -732,10 +954,23 @@ function l2_projection(q::TensorProductQuadrature{T, D}, f) where {T, D} end function l2_projection!(û::AbstractArray, q::TensorProductQuadrature{T, D}, f) where {T, D} - size(û) == size(q.basis) || throw(DimensionMismatch( - "the coefficient array is $(size(û)) but the basis is $(size(q.basis))")) + # 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) - copyto!(û, contract(q, F)) - mass_solve!(û, mass_operator(q), û) + l2_projection!(û, q, F, contract_workspace(q, eltype(û))) +end + +function l2_projection!(û::AbstractArray, q::TensorProductQuadrature{T, D}, + 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(q::TensorProductQuadrature, û::AbstractArray) + 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/Project.toml b/test/Project.toml index 143e2b3..4335039 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,6 +2,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FFTA = "b86e33f2-c0db-4aa1-a6e0-ab43e668529e" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -15,6 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Aqua = "0.8" Documenter = "1" FFTA = "0.3.2" +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..bcec4c9 --- /dev/null +++ b/test/integration/workspaces.jl @@ -0,0 +1,684 @@ +using ForwardDiff: 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` + @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 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 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") + # 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 + +@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 + +@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("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) + + û = 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))") + # 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 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 + +@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))) 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 + # (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 `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 = 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 = 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 = 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)") + # 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))") + # 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 + +# 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))") + # 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 + 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)) 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 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) + @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 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 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")