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..3d04374 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -157,11 +157,15 @@ 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} M::MT fact::FT + buf::Vector{T} function BandedMass(M::MT) where {T, MT <: AbstractMatrix{T}} kd = _bandwidth(M) @@ -170,7 +174,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 +304,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,10 +319,11 @@ 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 - # alternative is a plan that works everywhere except where it is passed a view. + # 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 -- @@ -329,7 +335,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 +480,41 @@ function mass_solve!(y::AbstractVector, op::FactorizedMass{T, MT, FT, :project}, return y end -# 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) +# 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. 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 + op.buf .= x + y .= ldiv!(op.fact, op.buf) + 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 + op.rbuf .= x + mul!(op.buf, op.plan, op.rbuf) + end op.buf .*= op.ĉ⁻¹ - mul!(y, op.iplan, op.buf) + if _contiguous(y) + mul!(y, op.iplan, op.buf) + else + 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..bdd99cd 100644 --- a/test/mass.jl +++ b/test/mass.jl @@ -312,6 +312,76 @@ 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 + + # 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) + 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 + @testset "$(rpad("accessors and errors",76))" begin q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16, 2π), 3)) op = mass_operator(q) 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