[michakraus]

Closes #10.

A caller that can describe a circulant matrix in O(n) no longer has to materialise O(n²)
entries to get its mass operator. The four points of the issue, additively — no existing
signature and no existing result changes.

What it adds

SimpleSplines.Circulant{T, CT}	an AbstractMatrix holding only the first column: C[k, l] = c[mod1(k - l + 1, n)], plus size and the generic Matrix, which materialises on demand
_check_circulant(::Circulant, …)	returns nothing — circulance holds by construction, so the quadratic verification has nothing left to verify. This is what keeps the path O(n) in construction as well as in storage
CirculantMass(c::AbstractVector, n)	wraps c and runs the existing constructor. mass_matrix(op) returns the Circulant; size, Matrix and every solve are unchanged
mass_operator(c::AbstractVector, b::PeriodicBSplineBasis)	the entry point, uniform mesh only

The case it comes from

The periodic B-spline Poisson solve in PoissonSolvers.jl shifts the periodic stiffness matrix
by the rank-one mean projector before factorising. S is sparse and banded, S + 𝟙𝟙ᵀ/n is
structurally full, and both are circulant.

At order 5 on a UniformMesh of 1024 cells, measured on Julia 1.13.0:

what the operator is handed	stored entries	summarysize
S — for reference	9 216	155 816 B
S .+ inv(n) — a SparseMatrixCSC, structurally full	1 048 576	16 785 576 B
S[:, 1] .+ inv(n) — a SparseVector	1 024	16 488 B
Vector(S[:, 1]) .+ inv(n) — a dense column	1 024	8 240 B

Three orders of magnitude, and the solves are identical to the last bit. The caller now writes
mass_operator(S[:, 1] .+ inv(n), b). Circulant is parametrised on the vector type, so a
sparse column is stored as one rather than being copied; Vector(…) first is half the size
again, if the caller cares.

What raises

Only a periodic basis on a UniformMesh has a circulant mass matrix, so that is the only
combination the vector form accepts. A periodic basis on a GradedMesh or a RandomMesh
raises an ArgumentError rather than falling to the FactorizedMass branch, which would
factorise a matrix the caller never described. A bounded basis matches no method at all.

Two decisions worth a reviewer's eye

• Circulant is not exported. The issue calls it an internal representation, and
  ToeplitzMatrices.jl exports the same name — using SimpleSplines, ToeplitzMatrices would
  then be an ambiguity error for a type nobody needs to name. It is documented on the Library
  page, and its docstring says how to reach it. Say the word if you would rather export it.
• CirculantMass(c, n) takes no atol. There is nothing on this path for a tolerance to
  govern, so the keyword is absent rather than accepted and ignored.

Verification

• Pkg.test() on this branch: 64 206 passed, 0 failed, Aqua.test_all included
  (ambiguities, piracy, stale deps, persistent tasks, compat bounds all green).
• Docs build (julia --project=docs -e 'include("docs/make.jl")'): exit 0, doctests ran,
  no missing-docstring and no broken-@ref warnings.
• JuliaFormatter --check, fatou lint and the NFC check: clean on every changed file.

julia-branch-verifier ran the seven-point pass against the diff. What it found:

Fixed

• test/mass_tests.jl:93 — a comment narrated the change instead of the state
  (circulance is no longer checked). Rewritten in the present tense.
• src/mass.jl:189 — # Verify rather than assume. was made incomplete by this branch, since
  one constructor path now assumes, justified by the type. The comment says so.
• src/mass.jl:216 — atol was accepted and ignored on the vector form. Removed, so it now
  raises a MethodError.

Raised and not accepted

• It reported the 16.8 MB figure in the changelog as unreproducible, measuring 8.39 MB. It
  had measured Matrix(S) .+ inv(n) — a dense Matrix{Float64}, 8 388 656 B. The expression
  the changelog names is S .+ inv(n), which stays sparse and so carries rowval beside
  nzval: SparseMatrixCSC, 1 048 576 stored entries, 16 785 576 B. The figure stands, and
  the three others in that sentence reproduce exactly.

Checked and clean

• Correctness: the mod1 indexing agrees with the documented C_{kl} formula and
  C[:, 1] == c. The rfft solve is correct for a general real circulant, not only a
  symmetric one — checked against a non-symmetric case to 1.4e-17. Out-of-range indices in
  all five corners raise BoundsError; linear and cartesian indexing agree.
• Consumers of the changed mass_matrix return type: *, mul! (including into a view),
  kron, copy and a two-factor KroneckerMass all work with a Circulant factor —
  km \ (km * B) round-trips to 6.4e-15.
• Type piracy: every new method is on an owned type or an owned name; no signature changed.
• Duplicated dependency functionality: no dependency of this package provides a circulant
  matrix type.
• Type stability: code_typed on every changed symbol returns concrete types; no Any
  originates in a changed hunk.
• Allocations, fresh process, --check-bounds=auto: a full getindex loop allocates 0 bytes;
  the vector construction path costs 32 B more than the matrix path at n = 1024 — the
  wrapper — the rest being the FFTW plans both paths share.

Not fixed

CirculantMass built from a materialised matrix still stores the whole matrix. Dropping M
there would make every periodic operator O(n), but it changes what mass_matrix returns for
existing callers, so it belongs in a separate decision, as the issue says.

🤖 Generated with Claude Code

[Copilot]

🟡 Changes recommended

Circulant mishandles valid offset-indexed vectors during indexing or materialization.

Get a fresh assessment by requesting another Copilot review.

Pull request overview

Adds an O(n)-storage and construction path for circulant mass operators from their first column.

Changes:

• Introduces internal Circulant storage and vector-based constructors.
• Adds uniform periodic mass_operator dispatch.
• Updates tests, documentation, and changelog coverage.

File summaries

File	Summary
test/mass_tests.jl	Tests construction, solving, storage, and rejection cases.
src/mass.jl	Implements circulant storage and construction; offset-indexed vectors require correction.
docs/src/usage/quadrature.md	Documents first-column usage.
docs/src/library.md	Documents Circulant.
CHANGELOG.md	Records the new API and limitation.

Review details

• Files reviewed: 5/5 changed files
• Comments generated: 1
• Review effort level: Lite

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[michakraus]

Closing this unmerged, in favour of #12, which removes the problem rather than making it cheaper
to carry.

Why

This branch exists because a caller had to materialise O(N²) entries to build the operator of a
circulant matrix it could describe in O(N). The caller in question is the periodic B-spline
Poisson solve, and what it was materialising was S + 𝟙𝟙ᵀ/N — the rank-one shift that makes the
singular periodic stiffness matrix invertible.

S + 𝟙𝟙ᵀ/N is circulant like S itself, so its whole effect is visible in the transform. It
changes one Fourier coefficient, from zero to one:

n= 32  ‖ℱ(c + 𝟙/N) - ℱ(c) - e₀‖∞ = 2.78e-17
n=128  ‖ℱ(c + 𝟙/N) - ℱ(c) - e₀‖∞ = 5.20e-18
n=512  ‖ℱ(c + 𝟙/N) - ℱ(c) - e₀‖∞ = 3.47e-18

The O(N²) fill-in buys exactly that one number. #12 sets the same coefficient without forming
anything, by deflating the kernel in the solve — kernel = :project, which gives the constant
mode a zero factor instead of an infinite one. That is the rule PoissonSolverFFT already applies
when it sets its k = 0 factor to zero.

With that in place the caller has no shift to build. It hands over the sparse assembly itself, and
nothing on the path is O(N²) any more: JuliaPlasma/PoissonSolvers.jl#10 takes a whole solver at
degree 4 on 2048 cells from 72 551 536 B to 5 426 160 B.

What that leaves

The first-column entry point is still a fair API for a caller that describes a circulant matrix in
O(N) and is not coming from a singular assembly. No such caller exists here now. What it would
save is the matrix a CirculantMass stores, which is the sparse banded assembly: 155 816 B at
degree 4 on 1024 cells, against 2 714 608 B for the whole solver. A constant factor on about 5 %
of the object, rather than an order.

#10 stays open for its own deferred item — making CirculantMass store the column in every
case, so that what the operator holds does not depend on how the matrix arrived. That is a smaller
change than this one and is worth judging on its own merits.

For anyone who revives this

Two things the review turned up, both worth carrying over:

• Circulant defines no mul!. mass_matrix(op) * x therefore runs at O(N²) through
  generic getindex. Every consumer works, but a KroneckerMass factor would be slow.
• The name collides with ToeplitzMatrices.Circulant, which is why this branch did not export
  it.

The branch is kept. Nothing here was wrong — it answered the question as asked, and the
measurements in it stand.