mass_operator takes only a materialised matrix, so a circulant operator costs O(n²) to build

[michakraus]

The problem

mass_operator(M, b) accepts an AbstractMatrix. On a periodic basis over a UniformMesh it
builds a CirculantMass, which reads one column of M and works from its transform
afterwards. The rest of M is used once, by _check_circulant, and then only ever reached again
through mass_matrix, size and Matrix.

A caller that wants the operator of a circulant matrix it can describe in O(n) therefore has to
materialise O(n²) entries first, and CirculantMass keeps them for the life of the operator.

The case this comes from is a periodic B-spline Poisson solve in
PoissonSolvers.jl. Constants lie in the kernel
of the periodic stiffness matrix, so the solver shifts it by the rank-one mean projector before
factorising:

S = stiffness_matrix(SplineQuadrature(b))   # sparse, banded, O(n) nnz
mass_operator(S .+ inv(size(S, 1)), b)      # structurally full, O(n²) nnz

S + 𝟙𝟙ᵀ/n is circulant, like S itself, and its first column is S[:, 1] .+ inv(n). The
solver knows that. It has no way to say it.

Measurements

The rank-one shift of the periodic stiffness matrix, order 5, on a UniformMesh of n cells.
Julia 1.13.0, SimpleSplines 0.1.0, macOS aarch64.

n	nnz(S)	nnz(S .+ inv(n))	summarysize
64	576	4 096	9.9 kB → 66 kB
1024	9 216	1 048 576	156 kB → 16.8 MB

The growth is O(n) to O(n²), so it is the size of the problem that decides how bad it gets. A
whole PoissonSolverSpline at order 5 on 2048 cells measures Base.summarysize 72 551 536 B,
and the CirculantMass reads 2048 numbers out of it.

The same solver against a Circulant of one column would be O(n) throughout.

Why the caller cannot work around it

• CirculantMass(M, n) is typed MT <: AbstractMatrix{T}, so a first column is not accepted.
• Writing a private circulant AbstractMatrix in the caller works for M[:, 1], but
  _check_circulant then falls to the generic method, which probes all n² positions one
  getindex at a time — the quadratic cost moves rather than going away.
• FactorizedMass avoids the fill-in only by giving up the FFT path.

Suggested fix

A Circulant matrix type carrying the first column, and one mass_operator method that takes a
vector.

1. Circulant{T} <: AbstractMatrix{T}, holding c::Vector{T}, with
   getindex(C, i, j) = C.c[mod1(i - j + 1, length(C.c))], size, and a Matrix that
   materialises on demand. About 15 lines. It is an internal representation, so it needs no
   arithmetic beyond indexing.

2. _check_circulant(::Circulant, …) = nothing. The property holds by construction, so the
   check has nothing left to verify. This is what makes the path O(n) end to end rather than
   O(n) in storage and O(n²) in construction.

3. CirculantMass(c::AbstractVector, n::Integer), which wraps c in a Circulant and runs
   the existing constructor. The stored M is then the Circulant, so mass_matrix(op) keeps
   returning an AbstractMatrix, size(op) keeps working, and Matrix(op) materialises only
   when a caller asks for it.

4. mass_operator(c::AbstractVector, b::PeriodicBSplineBasis), documented as taking the first
   column of a circulant matrix. It needs the UniformMesh branch only: a vector cannot be handed
   to FactorizedMass meaningfully, and a non-uniform mesh is not circulant, so that combination
   should raise rather than silently transpose into the other branch.

Nothing above changes an existing signature or an existing result, so it is additive.

What it does not fix

CirculantMass built from a materialised matrix still stores the whole matrix. Dropping M there
and reconstructing a Circulant from c would make every periodic operator O(n), but it changes
what mass_matrix returns for existing callers, so it belongs in a separate decision rather than
in this one.

Willing to open the PR

Happy to implement the four points above if the design suits you — say the word, or say which part
you would rather shape differently.