Solve a mass system on a non-contiguous view correctly

[michakraus]

Problem

mass_solve! on a BandedMass (every bounded basis) returns wrong values, with no error, when the result is a non-contiguous view such as view(w, 1:2:2m). It also writes into the parent entries that the view skips. The error depends on what the skipped entries hold. With 34 cells, a Dirichlet cubic basis and the skipped entries set to 7.0, the maximum error against a dense solve was about 2.8e3. ldiv!(y, op, x) and l2_projection! reach the same path.

Cause. BandedMass solves with ldiv! on a BandedMatrices banded Cholesky factor. For a strided vector that dispatches to pbtrs!, which calls LAPACK ?pbtrs with pointer(B) and does not check the stride of B (it calls chkstride1 on A only). LAPACK then reads and writes N contiguous entries from the first element of the view. A negative stride gives invalid argument #8 to LAPACK call.

The sibling CirculantMass fails loudly instead: an FFTW plan is bound to the strides it was planned for, so a strided result or right-hand side raises FFTW plan applied to wrong-strides array. FactorizedMass (CHOLMOD) already handled any stride, because it solves into a temporary and copies. KroneckerMass copies each fibre into contiguous buffers before the one-dimensional solve and was not affected.

Fix

BandedMass and CirculantMass hold a preallocated real buffer of length N. A non-contiguous argument (stride != 1, or not a StridedVector) goes through that buffer. A contiguous one takes the same path as before. Both solves stay allocation-free for any stride, negative strides included. The staging copies are broadcasts, so a strided argument of the wrong length raises DimensionMismatch, as the contiguous path does.

A copy is the right fix here, not a mask: both LAPACK ?pbtrs and an FFTW plan can only address memory with the layout they were given. The BLAS triangular solves accept a stride, but for a negative stride the pointer Julia passes is the wrong end of the array, and a probe segfaulted.

Test

New testset a non-contiguous argument is solved, not misread in test/mass.jl. It covers BandedMass (uniform Dirichlet, graded free), CirculantMass, FactorizedMass, and both periodic kernel = :project operators. Each case checks a stride-2 result, a stride-2 right-hand side, an aliased stride-2 argument and a negative-stride result against a dense solve. It also checks that the skipped parent entries stay unchanged, and that the banded and circulant solves allocate nothing on the strided path.

• On main: 6 pass, 8 fail, 7 error in that testset.
• On this branch: test/mass.jl gives 1089/1089, and the full Pkg.test() passes locally (Julia 1.13.1).

test/quality/jet.jl analyses mass_solve! at the stride-2 result type on both operators: 0 reports.

Pre-PR verification

Fixed

• The non-contiguous paths staged data with copyto!. That accepted a strided argument of the wrong length without an error, where the contiguous path throws. The staging copies are broadcasts, which raise DimensionMismatch. Two @test_throws cover it.
• A comment in the CirculantMass constructor is reflowed to the margin.
• test/quality/jet.jl analyses mass_solve! at the strided result type that the new @allocated test passes. The BandedMass docstring states that a non-contiguous argument goes through a buffer the operator owns.

Pre-existing, not changed here

• BandedMass on a contiguous result accepts a right-hand side that is too short (y === x || copyto!(y, x) in mass_solve!).
• The mass_solve! docstring says "two representations" where there are three.
• The missing stride check in BandedMatrices pbtrs! is an upstream defect; this PR does not depend on a fix there.

Checked and clean

• code_typed on mass_solve! for both operators at 5 argument shapes each: concrete return types, no dynamic dispatch.
• Zero allocations on six argument shapes per operator, fresh process, Julia 1.13.1.
• Aqua test_piracies, fatou lint and JuliaFormatter are clean. The CHANGELOG entry is present.

Not checked

• The Julia 1.11 floor locally (CI min covers it), a docs build, Float32 operators on the strided path, and thread safety.

🤖 Generated with Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[codecov]

Codecov Report

❌ Patch coverage is 94.11765% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 91.41%. Comparing base (e6999c4) to head (8076939).

Files with missing lines	Patch %	Lines
src/mass.jl	94.11%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #32      +/-   ##
==========================================
- Coverage   91.42%   91.41%   -0.02%     
==========================================
  Files           9        9              
  Lines        1178     1188      +10     
==========================================
+ Hits         1077     1086       +9     
- Misses        101      102       +1     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[michakraus]

JuliaDEC/SimpleSplines.jl PR #32: Solve a mass system on a non-contiguous view correctly

Verdict: request changes. (Posted as a comment review. The verdict is "request changes".)

The fix for the real-valued case is correct. I reproduced the defect and the fix. But the new BandedMass branch breaks a call that gave the correct answer before this PR: a complex non-contiguous result on a real operator now throws InexactError.

Findings

#	severity	file:line	claim	evidence
1	bug	src/mass.jl:493-499	Regression. A non-contiguous result whose element type is not the operator's T now goes through op.buf::Vector{T}. A ComplexF64 stride-2 result on a real BandedMass worked before this PR and now throws. The cause: ldiv!(op.fact, y) reaches LAPACK pbtrs! only when eltype(y) == T. Any other element type takes the generic LinearAlgebra path, which handles a stride correctly. Through l2_projection! this also breaks the complex samples that src/quadrature.jl:411-413 documents as accepted, when the result is a view on a bounded basis.	Probe at the PR head, Julia 1.13.1, a cubic Dirichlet basis with 34 cells. The pre-PR method body (copyto!(y, x); ldiv!(op.fact, y)) on a complex stride-2 y gives maxerr 4.564234742547372e-13, skipped entries changed = false. which returns ldiv!(C::Cholesky, B::AbstractVecOrMat) @ LinearAlgebra .../cholesky.jl:741. The PR head gives mass_solve! → throws InexactError: Float64(-1.4107270035180548 + 0.29255152095027787im). ldiv!(stride-2 complex y, opb, x) throws the same error. l2_projection!(stride-2 complex u, q, complex f) → throws InexactError: Float64(0.006251635923360947 - 0.012014747277459im).
2	quality	src/mass.jl:493	Proposed fix for #1. Use the buffer only where LAPACK is reached: if _contiguous(y) || eltype(y) !== T takes the direct path, and the else branch keeps the buffer (with op::BandedMass{T} and where {T} on the method).	I ran this body as a local function at the PR head. Stride-2, negative-stride and non-strided (view(_, collect(1:N))) results give maxerr 7.1e-14 for Float64 and maxerr 3.06e-13 for ComplexF64. The aliased results give the same values, and the skipped parent entries are untouched = true in both element types.
3	quality	test/mass.jl:315-375	The new testset uses only Float64 arguments, so it cannot see #1. Add a ComplexF64 stride-2 result on the banded cases, compared against the dense solve.	Read the testset: x = randn(N), and strided fills with 7.0.
4	quality	src/mass.jl:494	An adjacent pre-existing defect on a line that this PR moves. The PR body names it and leaves it. On the contiguous path, y === x || copyto!(y, x) accepts a right-hand side that is too short and solves on stale entries of y. The new strided branch on the next lines refuses the same input. The PR has made the two branches disagree on one input, so fix both in this change. A broadcast y .= x matches the new branch.	Probe: mass_solve!(zeros(N), op, ones(N - 1)) → returned, no error. mass_solve!(view(zeros(2N), 1:2:2N), op, ones(N - 1)) → throws DimensionMismatch.
5	nit	src/mass.jl:160-161	"so one operator is not to be shared between threads" is too broad for BandedMass. The contiguous branch reads only op.fact, and only a non-contiguous result writes op.buf. CirculantMass writes op.buf on every solve, but its docstring (src/mass.jl:269-271) has no thread note. Only the comment at src/tensorproduct.jl:413 says so.	Read src/mass.jl:492-517.
6	nit	src/mass.jl:450-451	The PR body names this one as pre-existing too: the mass_solve! docstring says "the two representations of MassOperator". There are four subtypes: FactorizedMass, BandedMass, CirculantMass and KroneckerMass. The docstring is in the function this PR changes.	grep -n '<: MassOperator' src/ returns src/mass.jl:67, :165, :301 and src/tensorproduct.jl:337.

Verified and fine

• I reproduced the defect on the pre-PR code paths, reached directly at the head. The banded stride-2 solve gives maxerr = 2769.3694192057733, skipped entries changed = true. The circulant stride-2 right-hand side gives ArgumentError: FFTW plan applied to wrong-strides array.
• Allocations, measured in fresh processes with Julia 1.13.1, --check-bounds=auto and --check-bounds=yes, minimum of 5 measurements. BandedMass and CirculantMass allocate 0 bytes each for these arguments: a contiguous one, a stride-2 result, a stride-2 right-hand side, an aliased stride-2 argument and a negative-stride result.
• Non-strided views (view(_, collect(1:N))) on both operators: maxerr 2.8e-13 and 3.8e-13 against the dense solve. op \ view(X, 1:2:2N, :) on BandedMass: maxerr 4.5e-13.
• A Float32 BandedMass with a Float64 result has the same accuracy on both branches: 7.0e-5 contiguous and 8.6e-5 strided. The CirculantMass element-type MethodErrors are on the contiguous path too, so this PR did not introduce them.
• The full Pkg.test() at the head passes on Julia 1.13.1. The "Mass operators" testset gives 1089/1089, and JET gives 12/12. opc in the new test/quality/jet.jl lines is defined at line 53.
• JuliaFormatter (sciml) reports the four changed .jl files as formatted.
• The PR body has no closing keyword, and closingIssuesReferences is empty.

CI

All required checks pass. Branch protection and the main ruleset require the same 7 contexts: Julia {min,1} - {ubuntu,macOS,windows}-latest - default and Doctests - ubuntu-latest. Neither the classic protection nor the ruleset requires Downgrade - ubuntu-latest. That job has continue-on-error: true (.github/workflows/CI.yml:128).

The Downgrade failure is a resolver floor check, not a test failure, and it is not caused by this PR. The PR changes no Project.toml. The job fails in the same way on main at e6999c4 (run 37038019344), which reports the workflow as success:

Error: forcedeps check failed: FFTW resolved to 1.3.1 but lower bound is 1.0.0
Error: forcedeps check failed: ContinuumArrays resolved to 0.20.5 but lower bound is 0.18.0
Error: forcedeps check failed: BandedMatrices resolved to 1.7.6 but lower bound is 1.0.0
ERROR: LoadError: forcedeps check failed for .: Some packages did not resolve to their lower bounds.

The suite never runs at the [compat] floors. This PR depends on BandedMatrices behaviour, so the floor BandedMatrices = "1" stays untested.

Not checked

• The PR body says the new testset gives "6 pass, 8 fail, 7 error" on main. I did not run the suite on main. I reproduced the defect on the pre-PR method bodies instead.
• Julia 1.11 locally. CI min is green.
• Thread safety beyond reading the code.

[michakraus]

Review findings: resolution

Follow-up to the review above. One commit pushed: 8076939 "Keep a non-T strided result of a BandedMass off the real buffer".

#	severity	finding	resolution
1	bug	A non-contiguous result of another element type than the factor's (ComplexF64 stride-2 on a real BandedMass) threw InexactError through op.buf::Vector{T}; main solved it correctly.	Fixed in 8076939. Reproduced first at f19b62f: InexactError(:Float64, ...); the pre-PR call ldiv!(op.fact, y) on the same input gave max error 4.6e-13 through ldiv!(C::Cholesky, B::AbstractVecOrMat) (LinearAlgebra cholesky.jl:741). After the fix: max error 3.3e-13 for stride 2, negative stride and aliased; skipped parent entries untouched; l2_projection! with complex samples into a stride-2 view matches the contiguous result exactly.
2	quality	Take the direct path when eltype(y) !== T.	Fixed in 8076939, as proposed. mass_solve!(y, op::BandedMass{T}, x) where {T} takes the direct path when _contiguous(y) || eltype(y) !== T. This matches the dispatch exactly: BandedMatrices sends only StridedVecOrMat{T} to pbtrs! (BandedCholesky.jl:90); every other element type goes to the generic solve. code_typed returns a concrete type for stride-2 Float64, stride-2 ComplexF64 and contiguous results.
3	quality	The new testset has no complex case.	Fixed in 8076939: a ComplexF64 stride-2 result on each BandedMass case, checked against the dense solve and the sentinel entries. test/mass.jl: 1093/1093 locally (Julia 1.13.1). The same call throws InexactError at f19b62f.
4	quality	The contiguous branch accepts a too-short right-hand side; the strided branch refuses it.	Not fixed. Pre-existing, and the PR body lists it as deliberately out of scope. One addition: the strided branch does not refuse a length-1 right-hand side either, because a broadcast expands it. mass_solve!(view(zeros(2N), 1:2:2N), op, [1.0]) and mass_solve!(zeros(N), op, [1.0]) both return with no error. An explicit length check on both branches would close both cases; that is a separate change.
5	nit	The thread note in the BandedMass docstring is broader than the code needs.	Not fixed. The note is conservative but correct: a solve on a non-contiguous argument writes op.buf, and a caller does not always know which path an argument takes. The missing note in the CirculantMass docstring is pre-existing and outside the changed hunks.
6	nit	The mass_solve! docstring says "two representations".	Not fixed. Pre-existing, outside the changed hunks, and named in the PR body.

CI at 8076939

• All seven required checks pass. Both classic branch protection and ruleset "main" require the same seven: Julia min and Julia 1 on ubuntu, macOS and windows, and Doctests - ubuntu-latest.
• Julia pre, Julia nightly, Documentation, codecov/patch and codecov/project also pass.
• Downgrade - ubuntu-latest fails, as it does on main (run 37038019344 at e6999c4). It is a [compat] floor check that fails before any test runs: forcedeps check failed: FFTW resolved to 1.3.1 but lower bound is 1.0.0, and likewise for ContinuumArrays (0.20.5 against 0.18.0) and BandedMatrices (1.7.6 against 1.0.0). This PR changes no Project.toml, and the job has continue-on-error.

No changelog change: the existing Unreleased entry already describes the final behaviour, and the regression in finding 1 was never released.