Build a CirculantMass from its first column

CHANGELOG.md

@@ -5,6 +5,37 @@ 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.2.0
+
+### Circulant mass matrices from first column alone
+
+A caller that can describe a circulant matrix in `O(n)` can now build its mass operator
+without materialising the `O(n²)` entries.
+
+`mass_operator(c::AbstractVector, b::PeriodicBSplineBasis)` takes that first column and
+returns a `CirculantMass` ready to solve, and `CirculantMass(c, n)` does the same one level
+down. The column is wrapped in a new internal type, `SimpleSplines.Circulant`, which stores
+nothing else. Circulance then holds by construction, so `_check_circulant` has nothing left
+to verify — which is what keeps the path `O(n)` in construction as well as in storage, rather
+than `O(n)` in storage and `O(n²)` in the check. `mass_matrix` returns the `Circulant`, and
+`size`, `Matrix` and every solve are unchanged.
+
+The motivating case is a periodic B-spline Poisson solve, which 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 the
+shift takes the matrix from 9 216 to 1 048 576 stored entries, and from 156 kB to 16.8 MB, and
+a `CirculantMass` built from it kept all of that for the life of the operator while reading
+1024 numbers out of it. The caller now writes `mass_operator(S[:, 1] .+ inv(n), b)`.
+
+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, and a bounded basis matches no method at all.
+
+**Known limitation:** a `CirculantMass` built from a materialised matrix still stores the whole
+matrix. Dropping it there would make every periodic operator `O(n)`, but it would change what
+`mass_matrix` returns for existing callers, so it is left for a separate decision.
+
 ## [0.1.0] — 2026-09-07
 
 The first release with any implementation in it. Before this the package was a

docs/src/library.md

@@ -117,6 +117,7 @@ l2_projection!
 ```@docs
 MassOperator
 CirculantMass
+Circulant
 BandedMass
 FactorizedMass
 KroneckerMass

docs/src/usage/quadrature.md

@@ -259,6 +259,38 @@ catch err
 end
 ```
 
+### When the caller already knows the matrix
+
+A caller that can describe a circulant matrix in ``O(N)`` need not materialise it. Hand
+[`mass_operator`](@ref) the **first column** instead of the matrix. The verification above then
+has nothing to check — circulance holds by construction — so storage and construction are both
+``O(N)`` rather than ``O(N^2)``.
+
+The case this exists for is a periodic Poisson solve. Constants lie in the kernel of the
+periodic stiffness matrix, so the solver shifts it by the rank-one mean projector before
+factorising. ``\mathbb{S}`` is sparse and banded, ``\mathbb{S} + \mathbb{1}\mathbb{1}^T/N`` is
+structurally full, and both are circulant.
+
+```@example asm
+bp = BSplineBasis(UniformMesh(64, 0 .. 2π), 3, Periodic())
+qp = SplineQuadrature(bp)
+Sp = stiffness_matrix(qp)
+N = nbasis(bp)
+
+heavy = mass_operator(Matrix(Sp) .+ inv(N), bp)     # N² entries, assembled and kept
+lean = mass_operator(Sp[:, 1] .+ inv(N), bp)        # N entries, and nothing else
+
+rhs = basis_integrals(qp)
+maximum(abs, lean \ rhs .- heavy \ rhs),
+Base.summarysize(mass_matrix(lean)), Base.summarysize(mass_matrix(heavy))
+```
+
+`mass_matrix(op)` returns the column wrapped in a `SimpleSplines.Circulant`, and `Matrix(op)`
+materialises the ``N^2`` entries only when something asks for them. Only a periodic basis on a
+[`UniformMesh`](@ref) has a circulant mass matrix, so this is the only combination the form
+accepts: a periodic basis on any other mesh raises rather than taking the
+[`FactorizedMass`](@ref) branch, and a bounded basis matches no method.
+
 ### Solving
 
 ```julia

src/mass.jl

@@ -107,8 +107,41 @@ function _bandwidth(M::AbstractMatrix)
     return kd
 end
 
+@doc raw"""
+    Circulant(c)
+
+The circulant matrix whose first column is `c`, stored as that column and nothing else.
+
+```math
+C_{kl} = c_{(k - l) \bmod n + 1} ,
+```
+
+so every column is the one to its left, rotated down by one. Storage and construction are
+both ``O(n)``, which is what lets [`mass_operator`](@ref) and [`CirculantMass`](@ref) be
+handed the first column of a circulant matrix instead of the ``n^2`` entries of its
+materialisation.
+
+This is an internal representation of a matrix the caller already knows to be circulant, so
+it carries no arithmetic beyond `getindex`, `size` and a `Matrix` that materialises the
+``n^2`` entries on demand. It is not exported; reach it as `SimpleSplines.Circulant`, or let
+`mass_operator(c, basis)` build it.
+"""
+struct Circulant{T, CT <: AbstractVector{T}} <: AbstractMatrix{T}
+    c::CT
+end
+
+Base.size(C::Circulant) = (length(C.c), length(C.c))
+
+# `mod1` wraps any index at all, so an out-of-range one would silently return an entry
+# rather than throw. The bounds check is what keeps that from happening.
+function Base.getindex(C::Circulant, k::Int, l::Int)
+    @boundscheck checkbounds(C, k, l)
+    return @inbounds C.c[mod1(k - l + 1, length(C.c))]
+end
+
 @doc raw"""
     CirculantMass(M, n)
+    CirculantMass(c, n)
 
 The mass matrix of a uniform periodic mesh, represented by the eigenvalues of its circulant
 structure and a pair of planned real transforms.
@@ -131,6 +164,13 @@ does not.
 The first column is read from the assembled matrix rather than recomputed, and the
 construction checks that the matrix really is circulant — a silent mismatch here would give
 wrong answers on every mesh that is uniform by accident rather than by construction.
+
+Given a vector instead, `c` is taken as that first column and wrapped in a
+[`Circulant`](@ref). Circulance then holds by construction, so the check has nothing left to
+verify and the whole path — storage as well as construction — is ``O(n)``. That is the form
+for a caller that can describe its matrix in ``O(n)``, such as a periodic stiffness matrix
+shifted by the rank-one mean projector; `mass_matrix` then returns the `Circulant`, and
+`Matrix` materialises only when something asks for it.
 """
 struct CirculantMass{T, MT, PT, IT} <: MassOperator{T}
     M::MT
@@ -146,9 +186,10 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T}
             "the mass matrix is $(size(M, 1))×$(size(M, 2)) but n = $(n)"))
         c = Vector{T}(M[:, 1])
 
-        # Verify rather than assume. A matrix that is banded but not circulant would still
-        # produce plausible numbers through the transform, and the failure would show up
-        # much later as a wrong conservation law.
+        # Verify rather than assume, wherever the argument does not already carry the
+        # guarantee. A matrix that is banded but not circulant would still produce plausible
+        # numbers through the transform, and the failure would show up much later as a wrong
+        # conservation law.
         _check_circulant(M, c, Int(n), atol)
 
         # The plans are made UNALIGNED so that they accept any strided argument -- a view
@@ -173,6 +214,10 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T}
     end
 end
 
+# No `atol`: there is nothing here for a tolerance to govern, and a keyword that is accepted
+# and then ignored is worse than one that is absent.
+CirculantMass(c::AbstractVector, n::Integer) = CirculantMass(Circulant(c), n)
+
 function _not_circulant(atol)
     ArgumentError(
         "the mass matrix is not circulant to within $(atol); a CirculantMass is only valid on " *
@@ -190,6 +235,12 @@ function _check_circulant(M::AbstractMatrix{T}, c::Vector{T}, n::Int, atol) wher
     return nothing
 end
 
+# A `Circulant` is circulant by construction, so the check has nothing to verify. Without
+# this method the generic one above would probe all n^2 positions one `getindex` at a time,
+# which is the cost the representation exists to avoid -- storage would be O(n) and
+# construction still O(n^2).
+_check_circulant(::Circulant{T}, ::Vector{T}, ::Int, ::Any) where {T} = nothing
+
 # The sparse check visits only the stored entries. Within a column the row indices are
 # distinct and r = mod1(i-j+1, n) is a bijection on 1:n, so counting the stored entries whose
 # predicted value is above the tolerance and comparing that count with the number of such
@@ -297,3 +348,27 @@ mass_operator(M::AbstractMatrix, ::AbstractBSplineBasis) = BandedMass(M)
 function mass_operator(M::AbstractMatrix, b::PeriodicBSplineBasis)
     mesh(b) isa UniformMesh ? CirculantMass(M, nbasis(b)) : FactorizedMass(M)
 end
+
+@doc raw"""
+    mass_operator(c::AbstractVector, b::PeriodicBSplineBasis)
+
+Build the [`CirculantMass`](@ref) of the circulant matrix whose **first column** is `c`,
+without materialising the other ``n^2 - n`` entries.
+
+This is the entry point for a caller that can describe its matrix in ``O(n)`` — a periodic
+stiffness matrix shifted by the rank-one mean projector, say, whose first column is
+`S[:, 1] .+ inv(n)` and which is circulant like `S` itself. Handing that matrix over
+assembled costs ``O(n^2)`` to build and ``O(n^2)`` to keep; handing over the column costs
+neither. `c` is wrapped in a [`Circulant`](@ref), which is what the operator stores.
+
+Only a periodic basis on a [`UniformMesh`](@ref) has a circulant mass matrix. A periodic
+basis on any other mesh raises rather than falling to the [`FactorizedMass`](@ref) branch,
+which would factorise a matrix the caller never described, and a bounded basis matches no
+method at all — a first column says nothing about a matrix that is not circulant.
+"""
+function mass_operator(c::AbstractVector, b::PeriodicBSplineBasis)
+    mesh(b) isa UniformMesh || throw(ArgumentError(
+        "a first column describes a circulant matrix, but the mass matrix of a periodic " *
+        "basis on a $(nameof(typeof(mesh(b)))) is not circulant; pass the assembled matrix"))
+    return CirculantMass(c, nbasis(b))
+end

test/mass_tests.jl

@@ -57,6 +57,63 @@ using Test
         end
     end
 
+    @testset "$(rpad("a circulant operator built from its first column alone",76))" begin
+        for p in 1:4, n in (12, 15, 16, 32)
+
+            b = PeriodicBSplineBasis(UniformMesh(n, 2π), p)
+            M = Matrix(mass_matrix(SplineQuadrature(b)))
+            c = M[:, 1]
+
+            C = SimpleSplines.Circulant(c)
+            @test C isa AbstractMatrix{Float64}
+            @test size(C) == (n, n)
+            @test Matrix(C) ≈ M
+            @test_throws BoundsError C[n + 1, 1]
+
+            op = mass_operator(c, b)
+            @test op isa CirculantMass
+            @test mass_matrix(op) isa SimpleSplines.Circulant
+            @test size(op) == (n, n)
+            @test Matrix(op) ≈ M
+
+            # the column and the whole matrix describe the same operator
+            ref = CirculantMass(M, n)
+            x = randn(n)
+            @test op \ x ≈ ref \ x atol = 1e-10
+            @test M * (op \ x) ≈ x atol = 1e-10
+        end
+    end
+
+    @testset "$(rpad("what the first-column form rejects",76))" begin
+        n = 16
+        b = PeriodicBSplineBasis(UniformMesh(n, 2π), 3)
+        c = Matrix(mass_matrix(SplineQuadrature(b)))[:, 1]
+
+        @test_throws DimensionMismatch CirculantMass(c, n - 1)
+        # this path does not check circulance, but it does check invertibility
+        @test_throws ArgumentError mass_operator(zeros(n), b)
+        # a periodic basis on a non-uniform mesh has no circulant mass matrix, so a first
+        # column describes nothing; it must not fall to the FactorizedMass branch instead
+        @test_throws ArgumentError mass_operator(c, PeriodicBSplineBasis(GradedMesh(n, 2π), 3))
+        @test_throws MethodError mass_operator(c, BSplineBasis(UniformMesh(n, 2π), 3))
+    end
+
+    @testset "$(rpad("the rank-one shift stays O(n) end to end",76))" begin
+        # the case this exists for: S is circulant and sparse, and S + 𝟙𝟙ᵀ/n is circulant
+        # and structurally full, so only its first column is worth carrying
+        n = 32
+        b = PeriodicBSplineBasis(UniformMesh(n, 2π), 3)
+        S = stiffness_matrix(SplineQuadrature(b))
+
+        dense = mass_operator(Matrix(S) .+ inv(n), b)
+        lean = mass_operator(S[:, 1] .+ inv(n), b)
+
+        x = randn(n)
+        @test lean \ x ≈ dense \ x atol = 1e-10
+        @test Matrix(lean) ≈ Matrix(dense)
+        @test Base.summarysize(mass_matrix(lean)) < Base.summarysize(mass_matrix(dense))
+    end
+
     @testset "$(rpad("in-place solve, aliasing, and no allocation",76))" begin
         q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(32, 2π), 3))
         op = mass_operator(q)