Solve a mass system on a non-contiguous view correctly

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

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
 

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!`

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)

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