From 1217223b2d61bb9f5a11fcbcab5b22546b25cc8c Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sat, 3 Oct 2026 20:24:09 +0200 Subject: [PATCH 1/4] Solve a mass system on a non-contiguous view correctly The LAPACK banded Cholesky solve behind BandedMass addresses its argument as contiguous memory, and the BandedMatrices wrapper does not check the stride, so mass_solve! into a stride-2 view returned wrong values and wrote into the parent entries the view skips. CirculantMass refused such a view through FFTW's stride check. Both operators now stage a non-contiguous argument through a buffer they own, which keeps the solve allocation-free. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 14 +++++++++++ src/mass.jl | 39 +++++++++++++++++++++++------- src/tensorproduct.jl | 11 ++++----- test/mass.jl | 57 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 107 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7c28a20..27d2a8c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,20 @@ 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.2 + +### Bug Fixes + +`mass_solve!`, `ldiv!` and `l2_projection!` on a bounded basis returned wrong values, with no +error, when the result was a non-contiguous view such as `view(w, 1:2:2m)`. They also wrote +into the entries of the parent array that the view skips. The cause is the LAPACK banded +Cholesky solve behind `BandedMass`: it addresses its argument as contiguous memory, and the +BandedMatrices wrapper does not check the stride. A `CirculantMass` refused such a view, or a +strided right-hand side, with FFTW's "plan applied to wrong-strides array". Both operators now +stage a non-contiguous argument through a buffer they own. The solve is correct for any stride, +negative strides included, and stays allocation-free. A contiguous argument takes the same path +as before. + ## [0.3.1] — 2026-10-02 ### New Features diff --git a/src/mass.jl b/src/mass.jl index 0068bb3..a71780e 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -162,6 +162,7 @@ the uniform periodic one. struct BandedMass{T, MT, FT} <: MassOperator{T} M::MT fact::FT + buf::Vector{T} function BandedMass(M::MT) where {T, MT <: AbstractMatrix{T}} kd = _bandwidth(M) @@ -170,7 +171,7 @@ struct BandedMass{T, MT, FT} <: MassOperator{T} issuccess(F) || throw(ArgumentError( "the mass matrix is not positive definite; the quadrature is too coarse to " * "resolve the basis")) - new{T, MT, typeof(F)}(M, F) + new{T, MT, typeof(F)}(M, F, Vector{T}(undef, size(M, 1))) end end @@ -300,6 +301,7 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} plan::PT iplan::IT buf::Vector{Complex{T}} + rbuf::Vector{T} n::Int function CirculantMass(M::MT, n::Integer; rtol = sqrt(eps(T)), @@ -314,9 +316,9 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} # much later as a wrong conservation law. _check_circulant(M, c, Int(n), rtol) - # The plans are made UNALIGNED so that they accept any strided argument -- a view + # The plans are made UNALIGNED so that they accept any contiguous argument -- a view # into a column of a matrix, in particular, whose alignment an aligned plan would - # reject at run time. At these sizes the difference is not measurable, and the + # reject at run time. Other strides go through `rbuf`; see `mass_solve!`. At these sizes the difference is not measurable, and the # alternative is a plan that works everywhere except where it is passed a view. # # ESTIMATE has to be given explicitly alongside it. Passing UNALIGNED alone replaces @@ -329,7 +331,8 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} iplan = plan_irfft(buf, n; flags = FFTW.ESTIMATE | FFTW.UNALIGNED) ĉ⁻¹ = _reciprocal_eigenvalues(plan * c, Int(n), kernel) - new{T, MT, typeof(plan), typeof(iplan)}(M, ĉ⁻¹, plan, iplan, buf, Int(n)) + new{T, MT, typeof(plan), typeof(iplan)}(M, ĉ⁻¹, plan, iplan, buf, + Vector{T}(undef, n), Int(n)) end end @@ -473,17 +476,37 @@ function mass_solve!(y::AbstractVector, op::FactorizedMass{T, MT, FT, :project}, return y end +# LAPACK's banded solve and FFTW's plans address their argument as contiguous memory. FFTW +# refuses any other stride; the banded LAPACK wrapper does not check, so a stride-2 view +# returns wrong values and writes into the parent entries the view skips. A non-contiguous +# argument therefore goes through a buffer the operator owns, which keeps both solves +# allocation-free on it, and a contiguous one goes straight through. +_contiguous(v::StridedVector) = stride(v, 1) == 1 +_contiguous(::AbstractVector) = false + # The banded factor does implement an in-place `ldiv!`, so this path allocates nothing. function mass_solve!(y::AbstractVector, op::BandedMass, x::AbstractVector) - y === x || copyto!(y, x) - ldiv!(op.fact, y) + if _contiguous(y) + y === x || copyto!(y, x) + ldiv!(op.fact, y) + else + copyto!(y, ldiv!(op.fact, copyto!(op.buf, x))) + end return y end function mass_solve!(y::AbstractVector, op::CirculantMass, x::AbstractVector) - mul!(op.buf, op.plan, x) + if _contiguous(x) + mul!(op.buf, op.plan, x) + else + mul!(op.buf, op.plan, copyto!(op.rbuf, x)) + end op.buf .*= op.ĉ⁻¹ - mul!(y, op.iplan, op.buf) + if _contiguous(y) + mul!(y, op.iplan, op.buf) + else + copyto!(y, mul!(op.rbuf, op.iplan, op.buf)) + end return y end diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index b45a642..b49565d 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -403,12 +403,11 @@ end # is viewed as (before, N_k, after) and each fibre along the middle index is solved on its own. # # 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 encodes the strides -# of the array it was planned for, not merely its alignment -- so the plan behind a -# `CirculantMass` rejects such a fibre outright with "plan applied to wrong-strides array". -# Creating the plans UNALIGNED, which is what lets them accept a *contiguous* column view, -# does not help here. The two buffers 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. +# 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!` diff --git a/test/mass.jl b/test/mass.jl index c9033a3..fef3b7d 100644 --- a/test/mass.jl +++ b/test/mass.jl @@ -312,6 +312,63 @@ Random.seed!(0x2f7a91c4) end end + @testset "$(rpad("a non-contiguous argument is solved, not misread",76))" begin + # LAPACK's banded solve and FFTW's plans address their argument as contiguous memory. + # A stride-2 view handed to either must still give the solution, and must leave the + # entries of the parent that the view skips alone. The baseline is a dense solve, not + # the contiguous path of the same operator. + n = 34 + cases = Any[] + for b in (BSplineBasis(UniformMesh(n, 0 .. 1), 3, Dirichlet()), + BSplineBasis(GradedMesh(n, 0 .. 1), 2), + PeriodicBSplineBasis(UniformMesh(n, 2π), 3), + PeriodicBSplineBasis(GradedMesh(n, 2π), 3)) + M = mass_matrix(SplineQuadrature(b)) + push!(cases, (mass_operator(M, b), x -> Matrix(M) \ x)) + end + for meshtype in (UniformMesh, GradedMesh) + b = PeriodicBSplineBasis(meshtype(n, 2π), 3) + S = stiffness_matrix(SplineQuadrature(b)) + shifted = Matrix(S) .+ inv(n) + push!(cases, (mass_operator(S, b; kernel = :project), + x -> shifted \ (x .- sum(x) / n))) + end + + sentinel = 7.0 + function strided(v) + w = fill(sentinel, 2length(v)) + w[1:2:end] .= v + return view(w, 1:2:length(w)) + end + for (op, solve) in cases + N = size(op, 1) + x = randn(N) + expected = solve(x) + + y = strided(zeros(N)) # a strided result + @test mass_solve!(y, op, x) ≈ expected atol = 1e-10 + @test all(==(sentinel), parent(y)[2:2:end]) + + y = zeros(N) # a strided right-hand side + @test mass_solve!(y, op, strided(x)) ≈ expected atol = 1e-10 + + z = strided(x) # both, aliased + @test mass_solve!(z, op, z) ≈ expected atol = 1e-10 + @test all(==(sentinel), parent(z)[2:2:end]) + + r = zeros(N) # a negative stride + @test mass_solve!(view(r, N:-1:1), op, x) ≈ expected atol = 1e-10 + + # the banded and circulant solves stage a strided argument through a buffer the + # operator owns, so they stay allocation-free on it + if !(op isa FactorizedMass) + y = strided(zeros(N)) + mass_solve!(y, op, x) + @test (@allocated mass_solve!(y, op, x)) == 0 + end + end + end + @testset "$(rpad("accessors and errors",76))" begin q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16, 2π), 3)) op = mass_operator(q) From cf2feb34ec26baebba1a7a5e55a125f65f9b5f48 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sat, 3 Oct 2026 20:33:23 +0200 Subject: [PATCH 2/4] Refuse a strided mass_solve! argument of the wrong length The buffered path for a non-contiguous argument copied through copyto!, which truncates silently: a strided rhs one entry short on a CirculantMass, or a strided result one entry long on either operator, returned a result where the contiguous path throws. The staging copies are broadcasts now, which raise DimensionMismatch. Two assertions in the non-contiguous testset cover it. A comment line in CirculantMass is reflowed to the margin, with its added sentence moved after the sentence whose referent it had displaced. Co-Authored-By: Claude Opus 5.5 --- src/mass.jl | 13 ++++++++----- test/mass.jl | 5 +++++ 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/mass.jl b/src/mass.jl index a71780e..febd0ac 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -318,8 +318,9 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} # The plans are made UNALIGNED so that they accept any contiguous argument -- a view # into a column of a matrix, in particular, whose alignment an aligned plan would - # reject at run time. Other strides go through `rbuf`; see `mass_solve!`. At these sizes the difference is not measurable, and the - # alternative is a plan that works everywhere except where it is passed a view. + # reject at run time. At these sizes the difference is not measurable, and the + # alternative is a plan that works everywhere except where it is passed a view. An + # argument of any other stride goes through `rbuf`; see `mass_solve!`. # # ESTIMATE has to be given explicitly alongside it. Passing UNALIGNED alone replaces # the flags rather than adding to them, and FFTW's default rigor then measures -- @@ -490,7 +491,8 @@ function mass_solve!(y::AbstractVector, op::BandedMass, x::AbstractVector) y === x || copyto!(y, x) ldiv!(op.fact, y) else - copyto!(y, ldiv!(op.fact, copyto!(op.buf, x))) + op.buf .= x + y .= ldiv!(op.fact, op.buf) end return y end @@ -499,13 +501,14 @@ function mass_solve!(y::AbstractVector, op::CirculantMass, x::AbstractVector) if _contiguous(x) mul!(op.buf, op.plan, x) else - mul!(op.buf, op.plan, copyto!(op.rbuf, x)) + op.rbuf .= x + mul!(op.buf, op.plan, op.rbuf) end op.buf .*= op.ĉ⁻¹ if _contiguous(y) mul!(y, op.iplan, op.buf) else - copyto!(y, mul!(op.rbuf, op.iplan, op.buf)) + y .= mul!(op.rbuf, op.iplan, op.buf) end return y end diff --git a/test/mass.jl b/test/mass.jl index fef3b7d..58431dc 100644 --- a/test/mass.jl +++ b/test/mass.jl @@ -365,6 +365,11 @@ Random.seed!(0x2f7a91c4) y = strided(zeros(N)) mass_solve!(y, op, x) @test (@allocated mass_solve!(y, op, x)) == 0 + + # a strided argument of the wrong length is refused, not truncated + @test_throws DimensionMismatch mass_solve!(strided(zeros(N + 1)), op, x) + @test_throws DimensionMismatch mass_solve!(strided(zeros(N)), op, + strided(x[1:(N - 1)])) end end end From f19b62f3b9aa780e6e4cfdd453e9f6a60f7c3cfe Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sat, 3 Oct 2026 20:35:51 +0200 Subject: [PATCH 3/4] Analyse the strided mass_solve! with JET; document the buffer test/quality/jet.jl analyses every function that a test checks with @allocated, at the argument types that test passes. The new testset checks mass_solve! into a stride-2 view, so JET now analyses that type on a banded and a circulant operator. The BandedMass docstring says that a non-contiguous argument goes through a buffer the operator owns. Co-Authored-By: Claude Opus 5.5 --- src/mass.jl | 3 +++ test/quality/jet.jl | 8 ++++++++ 2 files changed, 11 insertions(+) diff --git a/src/mass.jl b/src/mass.jl index febd0ac..2b4e28d 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -157,6 +157,9 @@ CHOLMOD, the factor answers `ldiv!` **in place**: a solve allocates nothing at a what makes [`l2_projection!`](@ref) allocation-free on every bounded basis rather than only on the uniform periodic one. +A non-contiguous argument, which LAPACK cannot address, is staged through a buffer the +operator owns, so one operator is not to be shared between threads. + [`FactorizedMass`](@ref) remains for the periodic non-uniform case, where the matrix wraps. """ struct BandedMass{T, MT, FT} <: MassOperator{T} diff --git a/test/quality/jet.jl b/test/quality/jet.jl index 5f9cd2e..02f00d0 100644 --- a/test/quality/jet.jl +++ b/test/quality/jet.jl @@ -62,6 +62,14 @@ if isdefined(JET, :JET_AVAILABLE) ? JET.JET_AVAILABLE : JET.JET_LOADABLE @test isempty(JET.get_reports(JET.report_opt(mass_solve!, (Vector{Float64}, T, Vector{Float64}); target_modules = (SimpleSplines,)))) end + + # test/mass.jl: mass_solve! into a stride-2 view, on a banded and a circulant operator + S2 = typeof(view(zeros(4), 1:2:4)) + opb = mass_operator(SplineQuadrature(BSplineBasis(UniformMesh(34, 0 .. 1), 3, Dirichlet()))) + @test isempty(JET.get_reports(JET.report_opt(mass_solve!, + (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,)))) else @test_skip "JET does not work on this Julia version" # aviatesk/JET.jl#681 end From 8076939830652e2a54b1f0ff4a53738a168bbac1 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sat, 3 Oct 2026 21:15:11 +0200 Subject: [PATCH 4/4] Keep a non-T strided result of a BandedMass off the real buffer A non-contiguous result of another element type than the factor's, such as a ComplexF64 stride-2 view on a real BandedMass, went through the Vector{T} buffer and threw InexactError. LAPACK is reached only for StridedVecOrMat{T}; any other element type takes LinearAlgebra's generic Cholesky ldiv!, which handles any stride, as it did before this branch. l2_projection! with complex samples into a view reaches the same path. A complex strided case is added to the testset. Co-Authored-By: Claude Opus 5.5 --- src/mass.jl | 8 +++++--- test/mass.jl | 8 ++++++++ 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/mass.jl b/src/mass.jl index 2b4e28d..3d04374 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -488,9 +488,11 @@ end _contiguous(v::StridedVector) = stride(v, 1) == 1 _contiguous(::AbstractVector) = false -# The banded factor does implement an in-place `ldiv!`, so this path allocates nothing. -function mass_solve!(y::AbstractVector, op::BandedMass, x::AbstractVector) - if _contiguous(y) +# The banded factor does implement an in-place `ldiv!`, so this path allocates nothing. Only a +# result of the factor's own element type reaches LAPACK; any other one takes the generic +# triangular solve, which handles any stride and which the real buffer could not hold. +function mass_solve!(y::AbstractVector, op::BandedMass{T}, x::AbstractVector) where {T} + if _contiguous(y) || eltype(y) !== T y === x || copyto!(y, x) ldiv!(op.fact, y) else diff --git a/test/mass.jl b/test/mass.jl index 58431dc..bdd99cd 100644 --- a/test/mass.jl +++ b/test/mass.jl @@ -359,6 +359,14 @@ Random.seed!(0x2f7a91c4) r = zeros(N) # a negative stride @test mass_solve!(view(r, N:-1:1), op, x) ≈ expected atol = 1e-10 + # a complex strided result of a real banded operator, which LAPACK never sees + if op isa BandedMass + xc = randn(ComplexF64, N) + yc = view(fill(complex(sentinel), 2N), 1:2:(2N)) + @test mass_solve!(yc, op, xc) ≈ solve(xc) atol = 1e-10 + @test all(==(sentinel), parent(yc)[2:2:end]) + end + # the banded and circulant solves stage a strided argument through a buffer the # operator owns, so they stay allocation-free on it if !(op isa FactorizedMass)