diff --git a/CHANGELOG.md b/CHANGELOG.md index 8247a86..c835fa0 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/src/library.md b/docs/src/library.md index eaa5364..9ecb67c 100644 --- a/docs/src/library.md +++ b/docs/src/library.md @@ -117,6 +117,7 @@ l2_projection! ```@docs MassOperator CirculantMass +Circulant BandedMass FactorizedMass KroneckerMass diff --git a/docs/src/usage/quadrature.md b/docs/src/usage/quadrature.md index 1ffbf52..029ef86 100644 --- a/docs/src/usage/quadrature.md +++ b/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 diff --git a/src/mass.jl b/src/mass.jl index ae81dd6..b0b312f 100644 --- a/src/mass.jl +++ b/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 diff --git a/test/mass_tests.jl b/test/mass_tests.jl index fca917e..72927ac 100644 --- a/test/mass_tests.jl +++ b/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)