From e52e954ad1e94ae8bfbc2b299cb66f129a0d9740 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Fri, 21 Aug 2026 19:18:55 +0900 Subject: [PATCH 1/4] Make l2_projection! and basis_integrals allocation-free MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `l2_projection!` was in-place only in its output: it built the `f ⊙ w` product and the load vector as fresh arrays on every call, 6.3 kB per projection at N = 128. The quadrature now holds an `f ⊙ w` buffer and the load vector is formed with `mul!` straight into `û`, which takes a uniform-mesh projection to zero allocations end to end -- the `CirculantMass` solve was already allocation-free. On a non-uniform mesh the CHOLMOD temporary remains, CHOLMOD having no in-place `ldiv!`. The shared buffer makes the method non-reentrant across threads, so the docstring says so and points at the allocating `l2_projection` as the way out. It also gains the sample-length check the allocating method already had; the buffer would otherwise have turned a wrong-length `f` into a broadcast error rather than a named one. `basis_integrals` recomputed `Φ₀ * w` on every call, 1.1 kB a time, although ∫ φ_i dx is as much a constant of the discretisation as the mass matrix. It is now assembled with the quadrature and returned by reference, like `mass_matrix`, and the docstring records that it must not be mutated. Both properties are pinned by tests, in the style of the existing `mass_solve!` allocation test. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 12 +++++++++++- src/quadrature.jl | 30 +++++++++++++++++++++++++++--- test/quadrature_tests.jl | 15 +++++++++++++++ 3 files changed, 53 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 33cf2f3..6de4833 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -86,7 +86,17 @@ was the difference between 6.0 ms and 0.4 ms for assembling one Hessian. The constant assemblies — `mass_matrix`, `stiffness_matrix`, `derivative_matrix` and any `mixed_matrix` — are memoised on first use. They do not depend on the field, but a downstream -time integrator asks for them inside every Newton iteration of every step. +time integrator asks for them inside every Newton iteration of every step. `basis_integrals` +is one of these constants too, and is now assembled with the quadrature and returned by +reference rather than recomputed per call. + +The paths a time integrator actually runs per step allocate nothing. Evaluation is +allocation-free at any degree and derivative order on any mesh; a `CirculantMass` solve is +allocation-free; and `l2_projection!` is now allocation-free end to end on a uniform mesh, +the `f ⊙ w` product going into a buffer held by the quadrature and the load vector being +formed with `mul!` straight into the output. That buffer makes `l2_projection!` +non-reentrant, which its docstring says. On a non-uniform mesh a CHOLMOD temporary remains, +CHOLMOD having no in-place `ldiv!`. ### Dependencies diff --git a/src/quadrature.jl b/src/quadrature.jl index b17427a..465b296 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -78,6 +78,8 @@ struct SplineQuadrature{T, BT <: PeriodicBSplineBasis{T}, MO <: MassOperator{T}} w::Vector{T} Φ::Vector{SparseMatrixCSC{T, Int}} mass::MO + integrals::Vector{T} + scratch::Vector{T} cache::Dict{Tuple{Int,Int}, SparseMatrixCSC{T, Int}} function SplineQuadrature(basis::BT; @@ -149,7 +151,13 @@ struct SplineQuadrature{T, BT <: PeriodicBSplineBasis{T}, MO <: MassOperator{T}} rethrow() end - new{T, BT, typeof(mass)}(basis, Int(nq), x, w, Φ, mass, + # ∫ φ_i dx is a constant of the discretisation like the mass matrix, so it is + # assembled here rather than recomputed on every call. `scratch` is the f ⊙ w + # buffer that keeps `l2_projection!` from allocating one per call. + integrals = Φ[1] * w + scratch = Vector{T}(undef, n * nq) + + new{T, BT, typeof(mass)}(basis, Int(nq), x, w, Φ, mass, integrals, scratch, Dict{Tuple{Int,Int}, SparseMatrixCSC{T, Int}}()) end end @@ -306,8 +314,11 @@ This is the gradient of the total mass ``C_0 = \int_\Omega u \, dx`` with respec degrees of freedom, and it equals ``\mathbb{M} \mathbf{1}`` because the basis is a partition of unity. It spans the kernel of the first discrete bracket, which is why the mass is a Casimir there. + +Assembled once when the [`SplineQuadrature`](@ref) is built and returned by reference, as +[`mass_matrix`](@ref) is — do not mutate the result. """ -basis_integrals(q::SplineQuadrature) = basis_values(q, 0) * q.w +basis_integrals(q::SplineQuadrature) = q.integrals @doc raw""" l2_projection(q::SplineQuadrature, f) @@ -339,11 +350,24 @@ end l2_projection!(û, q::SplineQuadrature, f) In-place [`l2_projection`](@ref), writing the coefficients into `û`. + +On a [`UniformMesh`](@ref) this allocates nothing: the ``f \odot w`` product goes into a +buffer held by the quadrature, the load vector is formed with `mul!` straight into `û`, and +the [`CirculantMass`](@ref) solve is itself allocation-free. On a non-uniform mesh the +CHOLMOD solve still allocates a temporary, as [`mass_solve!`](@ref) notes. + +That shared buffer is also what makes this method non-reentrant: two threads projecting +through the same [`SplineQuadrature`](@ref) at once would overwrite each other's `f ⊙ w`. +Give each thread its own quadrature, or use the allocating [`l2_projection`](@ref). """ function l2_projection!(û::AbstractVector, q::SplineQuadrature, f::AbstractVector) length(û) == nbasis(q) || throw(DimensionMismatch( "the coefficient vector has $(length(û)) entries but the basis has $(nbasis(q))")) - û .= basis_values(q, 0) * (q.w .* f) + length(f) == length(q.x) || throw(DimensionMismatch( + "the function was sampled at $(length(f)) points but the quadrature has " * + "$(length(q.x))")) + q.scratch .= q.w .* f + mul!(û, basis_values(q, 0), q.scratch) mass_solve!(û, q.mass, û) return û end diff --git a/test/quadrature_tests.jl b/test/quadrature_tests.jl index 7ef5adc..20fe773 100644 --- a/test/quadrature_tests.jl +++ b/test/quadrature_tests.jl @@ -59,6 +59,12 @@ using Test @test mass_matrix(q) * ones(16) ≈ Iv @test all(>(0), Iv) end + # a constant of the discretisation, assembled once and returned by reference rather + # than recomputed on every call + q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16), 3)) + @test basis_integrals(q) === basis_integrals(q) + basis_integrals(q) + @test (@allocated basis_integrals(q)) == 0 end @testset "$(rpad("S is already antisymmetric, given nq >= p",76))" begin @@ -156,6 +162,15 @@ using Test @test l2_projection!(ŵ, q, sin) ≈ û @test_throws DimensionMismatch l2_projection(q, [1.0, 2.0]) @test_throws DimensionMismatch l2_projection!(zeros(3), q, sin) + @test_throws DimensionMismatch l2_projection!(similar(û), q, [1.0, 2.0]) + + # on a uniform mesh the whole projection is allocation-free: the f ⊙ w buffer lives + # in the quadrature, the load vector is formed with mul! straight into û, and the + # CirculantMass solve allocates nothing + fq = sin.(quadrature_nodes(q)) + l2_projection!(ŵ, q, fq) + @test (@allocated l2_projection!(ŵ, q, fq)) == 0 + @test ŵ ≈ û end @testset "$(rpad("L2 projection converges at order p+1",76))" begin From d82b74a828ff253ea18f4b5d7abfd3a50ccd6f49 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Fri, 21 Aug 2026 19:19:24 +0900 Subject: [PATCH 2/4] Modernise the CI workflows and add AUTHORS.md Split the documentation build out of `CI.yml` into its own `Documenter.yml`, so that a docs failure and a test failure are separate signals and the docs job is not queued behind the test matrix. The README gains a badge for it. Bring the actions up to current versions -- `actions/checkout@v7`, `julia-actions/setup-julia@v3`, `julia-actions/cache@v3`, `codecov/codecov-action@v7`, and `RegisterAction` pinned so Dependabot can bump it, which `@latest` prevented. Dependabot now groups action bumps into one pull request. The test matrix names Julia versions by alias rather than by number: `min` resolves the lower bound of the `julia` compat entry, so the matrix tracks the declared support window instead of having to be edited alongside it; `lts` and `1` cover the long-term-support and current stable releases; `pre` and `nightly` run on Linux only under `continue-on-error`, an upcoming-release failure being information rather than a broken build. `arch` becomes `default`, which tests Julia natively on the ARM64 macOS runners instead of under Rosetta. Two workflow bugs fixed on the way: - `docs/Project.toml` pinned `CompactBasisFunctions` to an absolute path on a developer's machine through a `[sources]` entry. That path does not exist on a runner, so the documentation build could only ever have succeeded locally. 0.3.1 is registered, so the entry goes and the dependency resolves from the registry. - `CompatHelper.yml` invoked `julia` without installing it. The runner images no longer ship a Julia, so the workflow failed before CompatHelper started. Add `AUTHORS.md` and rename `LICENSE` to `LICENSE.md`, whose copyright line now names "The SimpleSplines Authors" and points at it -- otherwise AUTHORS.md's reference to what the LICENSE mentions would not be true. Add the `.githooks/pre-push` hook the README already told the reader to enable. `.githooks/` did not exist, so `git config core.hooksPath .githooks` was silently a no-op. It runs the test suite and refuses the push if it fails, with `SIMPLESPLINES_SKIP_TESTS=1` as the documented override. Co-Authored-By: Claude Opus 5 (1M context) --- .githooks/pre-push | 39 +++++++++++++++++ .github/dependabot.yml | 5 +++ .github/workflows/CI.yml | 68 +++++++++++++----------------- .github/workflows/CompatHelper.yml | 10 +++++ .github/workflows/Documenter.yml | 48 +++++++++++++++++++++ .github/workflows/Register.yml | 5 ++- .github/workflows/TagBot.yml | 1 + AUTHORS.md | 18 ++++++++ CHANGELOG.md | 31 ++++++++++++++ LICENSE => LICENSE.md | 2 +- README.md | 13 ++++++ docs/Project.toml | 4 +- 12 files changed, 201 insertions(+), 43 deletions(-) create mode 100755 .githooks/pre-push create mode 100644 .github/workflows/Documenter.yml create mode 100644 AUTHORS.md rename LICENSE => LICENSE.md (93%) diff --git a/.githooks/pre-push b/.githooks/pre-push new file mode 100755 index 0000000..b9fc049 --- /dev/null +++ b/.githooks/pre-push @@ -0,0 +1,39 @@ +#!/bin/sh +# +# Run the test suite before a push, and refuse the push if it fails. +# +# Enable with +# +# git config core.hooksPath .githooks +# +# and skip a single push with `SIMPLESPLINES_SKIP_TESTS=1 git push` — for a +# documentation-only change, say, or when the failure is already known and being +# pushed to a branch on purpose. + +set -eu + +if [ "${SIMPLESPLINES_SKIP_TESTS:-0}" != "0" ]; then + printf 'pre-push: SIMPLESPLINES_SKIP_TESTS is set, skipping the test suite.\n' >&2 + exit 0 +fi + +if ! command -v julia >/dev/null 2>&1; then + printf 'pre-push: julia is not on PATH, cannot run the test suite; refusing the push.\n' >&2 + printf ' Set SIMPLESPLINES_SKIP_TESTS=1 to push anyway.\n' >&2 + exit 1 +fi + +# The hook runs with the working directory at the top level of the work tree already, but +# `git push` from a subdirectory of a worktree is not worth relying on that for. +root=$(git rev-parse --show-toplevel) + +printf 'pre-push: running the SimpleSplines test suite...\n' >&2 + +if julia --project="$root" -e 'using Pkg; Pkg.test()'; then + printf 'pre-push: tests passed.\n' >&2 + exit 0 +fi + +printf 'pre-push: tests FAILED; refusing the push.\n' >&2 +printf ' Set SIMPLESPLINES_SKIP_TESTS=1 to push anyway.\n' >&2 +exit 1 diff --git a/.github/dependabot.yml b/.github/dependabot.yml index 700707c..bebf22d 100644 --- a/.github/dependabot.yml +++ b/.github/dependabot.yml @@ -5,3 +5,8 @@ updates: directory: "/" # Location of package manifests schedule: interval: "weekly" + # One pull request for all action bumps rather than one per action. + groups: + github-actions: + patterns: + - "*" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 923306b..3d10615 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -16,65 +16,55 @@ jobs: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} runs-on: ${{ matrix.os }} timeout-minutes: 60 + # A pre-release or nightly failure is information, not a broken build. + continue-on-error: ${{ matrix.experimental }} permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created actions: write contents: read strategy: fail-fast: false matrix: + # 'min' resolves the lower bound of the julia compat entry in Project.toml, so the + # matrix follows the declared support window instead of having to be edited alongside + # it; '1' is the current stable release. version: - - '1.10' - - '1.12' - - 'nightly' + - 'min' + - '1' os: - ubuntu-latest - macOS-latest - windows-latest arch: - - x64 + # 'default' is the runner's native architecture, which is aarch64 on macOS-latest; + # a hardcoded 'x64' would test Julia under Rosetta there. + - default + experimental: + - false + include: + - version: 'lts' + os: ubuntu-latest + arch: default + experimental: false + - version: 'pre' + os: ubuntu-latest + arch: default + experimental: true + - version: 'nightly' + os: ubuntu-latest + arch: default + experimental: true steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 + - uses: actions/checkout@v7 + - uses: julia-actions/setup-julia@v3 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 + - uses: codecov/codecov-action@v7 with: files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} fail_ci_if_error: false - docs: - name: Documentation - runs-on: ubuntu-latest - permissions: - actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created - contents: write - statuses: write - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: '1' - - uses: julia-actions/cache@v2 - - name: Configure doc environment - shell: julia --project=docs --color=yes {0} - run: | - using Pkg - Pkg.develop(PackageSpec(path=pwd())) - Pkg.instantiate() - - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-docdeploy@v1 - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} - - name: Run doctests - shell: julia --project=docs --color=yes {0} - run: | - using Documenter: DocMeta, doctest - using SimpleSplines - DocMeta.setdocmeta!(SimpleSplines, :DocTestSetup, :(using SimpleSplines); recursive=true) - doctest(SimpleSplines) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index cba9134..469681a 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -3,10 +3,20 @@ on: schedule: - cron: 0 0 * * * workflow_dispatch: +permissions: + contents: write + pull-requests: write jobs: CompatHelper: runs-on: ubuntu-latest + timeout-minutes: 30 steps: + # The runner images no longer ship a Julia, so one has to be installed rather than + # assumed; without this the `julia -e` steps below fail before CompatHelper starts. + - uses: julia-actions/setup-julia@v3 + with: + version: '1' + - uses: julia-actions/cache@v3 - name: Pkg.add("CompatHelper") run: julia -e 'using Pkg; Pkg.add("CompatHelper")' - name: CompatHelper.main() diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml new file mode 100644 index 0000000..6902c2d --- /dev/null +++ b/.github/workflows/Documenter.yml @@ -0,0 +1,48 @@ +name: Documentation +on: + push: + branches: + - main + tags: ['*'] + pull_request: + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + build: + name: Documentation + runs-on: ubuntu-latest + timeout-minutes: 60 + permissions: + actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created + contents: write + pull-requests: read + statuses: write + steps: + - uses: actions/checkout@v7 + - uses: julia-actions/setup-julia@v3 + with: + version: '1' + - uses: julia-actions/cache@v3 + - uses: julia-actions/julia-buildpkg@v1 + # docs/Project.toml carries a [sources] entry for the package itself, so instantiating + # the docs environment is all that is needed to resolve it against the checkout. + - name: Instantiate the docs environment + shell: julia --project=docs --color=yes {0} + run: | + using Pkg + Pkg.instantiate() + - uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + - name: Run doctests + shell: julia --project=docs --color=yes {0} + run: | + using Documenter: DocMeta, doctest + using SimpleSplines + DocMeta.setdocmeta!(SimpleSplines, :DocTestSetup, :(using SimpleSplines); recursive=true) + doctest(SimpleSplines) diff --git a/.github/workflows/Register.yml b/.github/workflows/Register.yml index 5b7cd3b..0367532 100644 --- a/.github/workflows/Register.yml +++ b/.github/workflows/Register.yml @@ -8,9 +8,10 @@ on: jobs: register: runs-on: ubuntu-latest + timeout-minutes: 30 permissions: - contents: write + contents: write steps: - - uses: julia-actions/RegisterAction@latest + - uses: julia-actions/RegisterAction@v0.3.2 with: token: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index 0cd3114..234915d 100644 --- a/.github/workflows/TagBot.yml +++ b/.github/workflows/TagBot.yml @@ -24,6 +24,7 @@ jobs: TagBot: if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' runs-on: ubuntu-latest + timeout-minutes: 30 steps: - uses: JuliaRegistries/TagBot@v1 with: diff --git a/AUTHORS.md b/AUTHORS.md new file mode 100644 index 0000000..3489aaf --- /dev/null +++ b/AUTHORS.md @@ -0,0 +1,18 @@ +# Authors + +SimpleSplines' development is coordinated by a group of *principal developers*, who are also +its main contributors and who can be contacted in case of questions about SimpleSplines. In +addition, there are *contributors* who have provided substantial additions or modifications. +Together, these two groups form "The SimpleSplines Authors" as mentioned in the +[LICENSE](LICENSE.md) file. + +## Principal Developers + +* [Michael Kraus](https://www.michael-kraus.org/), + Max Planck Institute for Plasma Physics, Garching, Germany + +## Contributors + +The following people contributed to SimpleSplines and are listed in alphabetical order: + +* Michael Kraus diff --git a/CHANGELOG.md b/CHANGELOG.md index 6de4833..fad089b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -104,9 +104,40 @@ CHOLMOD having no in-place `ldiv!`. assemblies do not map onto that package's Galerkin interface. `FFTW` and `SparseArrays` were added for the above, and `CompactBasisFunctions` for the shared `Basis` hierarchy. +### Repository and CI + +The documentation build moved out of `CI.yml` into its own `Documenter.yml` workflow, so that +a docs failure and a test failure are separate signals and the docs job is not queued behind +the test matrix. The README gained a badge for it. + +The workflows were brought up to current action versions — `actions/checkout@v7`, +`julia-actions/setup-julia@v3`, `julia-actions/cache@v3`, `codecov/codecov-action@v7` — and +the test matrix now names Julia versions by alias rather than by number: + +- `min` resolves the lower bound of the `julia` compat entry, so the matrix tracks the + declared support window instead of having to be edited alongside it; +- `lts` and `1` cover the long-term-support and current stable releases; +- `pre` and `nightly` run on Linux only and are `continue-on-error`, since an upcoming-release + failure is information rather than a broken build. + +`arch` is now `default` rather than `x64`, which is what tests Julia natively on the ARM64 +macOS runners instead of under Rosetta. + +`AUTHORS.md` was added and `LICENSE` renamed to `LICENSE.md`, whose copyright line now names +"The SimpleSplines Authors" and points at it. + +The `pre-push` hook the README asks the reader to enable now exists in `.githooks/`. It runs +the test suite and refuses the push if it fails; `SIMPLESPLINES_SKIP_TESTS=1` overrides it. + ### Fixed - `QuadratureRules` compat was `"0.1"`, which could not co-resolve with `CompactBasisFunctions`; it is now `"0.2"`. - `LinearAlgebra` compat was `"1.12.0"` alongside `julia = "1.10"`, which contradicted the 1.10 row of the CI matrix; it is now `"1"`. +- `docs/Project.toml` pinned `CompactBasisFunctions` to an absolute path on a developer's + machine through a `[sources]` entry. That path does not exist on a CI runner, so the + documentation build could only ever have succeeded locally; the entry is removed and the + dependency now resolves from the registry. +- `CompatHelper.yml` invoked `julia` without installing it. The runner images no longer ship + a Julia, so the workflow failed before CompatHelper started; it now sets Julia up first. diff --git a/LICENSE b/LICENSE.md similarity index 93% rename from LICENSE rename to LICENSE.md index b472609..7102389 100644 --- a/LICENSE +++ b/LICENSE.md @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2025 Michael Kraus +Copyright (c) The SimpleSplines Authors (see [AUTHORS.md](AUTHORS.md)) Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/README.md b/README.md index 0bae0c0..b8f18ed 100644 --- a/README.md +++ b/README.md @@ -3,6 +3,7 @@ [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaDEC.github.io/SimpleSplines.jl/stable/) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaDEC.github.io/SimpleSplines.jl/dev/) [![Build Status](https://github.com/JuliaDEC/SimpleSplines.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaDEC/SimpleSplines.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Documentation](https://github.com/JuliaDEC/SimpleSplines.jl/actions/workflows/Documenter.yml/badge.svg?branch=main)](https://github.com/JuliaDEC/SimpleSplines.jl/actions/workflows/Documenter.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/JuliaDEC/SimpleSplines.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaDEC/SimpleSplines.jl) [![PkgEval](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/S/SimpleSplines.svg)](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/S/SimpleSplines.html) @@ -19,3 +20,15 @@ To run the test suite before every push, enable the repository's git hooks: ```sh git config core.hooksPath .githooks ``` + +`.githooks/pre-push` refuses the push if the suite fails. To push regardless — a +documentation-only change, or a failure that is already known — set + +```sh +SIMPLESPLINES_SKIP_TESTS=1 git push +``` + +## License + +SimpleSplines is licensed under the [MIT License](LICENSE.md). See [AUTHORS.md](AUTHORS.md) +for the list of authors the copyright refers to. diff --git a/docs/Project.toml b/docs/Project.toml index ef2915a..7b31ff2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,5 +4,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" SimpleSplines = "ac480f6c-9a94-4d53-a283-9db69f109154" [sources] -CompactBasisFunctions = {path = "/Users/mkraus/Datashare/Julia/CompactBasisFunctions"} SimpleSplines = {path = ".."} + +[compat] +Documenter = "1" From 80866d0a8ffb87a2a39aaa92609d2dedccbf7a35 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Fri, 21 Aug 2026 19:22:23 +0900 Subject: [PATCH 3/4] Record the weighted_matrix allocation as an open issue `weighted_matrix` is the one assembly that cannot be memoised, since it depends on the field, so a time integrator asks for a different one inside every Newton iteration -- the call pattern under which allocation matters most. It costs 260 kB a call at N = 128, p = 3, and only 5 kB of that is a temporary; the rest is a freshly built SparseMatrixCSC whose structure is recomputed each time even though the sparsity pattern of Phi_a diag(f w) Phi_b' does not depend on f. Recorded rather than fixed: the fix needs a `weighted_matrix!` entry point, because the present signature has nowhere to write, and the sparse triple product has to be written out by hand against a cached pattern instead of delegating to SparseArrays. That widens the API, so it is the caller's call. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fad089b..726bef8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -141,3 +141,37 @@ the test suite and refuses the push if it fails; `SIMPLESPLINES_SKIP_TESTS=1` ov dependency now resolves from the registry. - `CompatHelper.yml` invoked `julia` without installing it. The runner images no longer ship a Julia, so the workflow failed before CompatHelper started; it now sets Julia up first. + +## Open Issues + +### `weighted_matrix` allocates a fresh matrix per call + +`weighted_matrix` is the one assembly that cannot be memoised — it depends on the field, so a +downstream time integrator asks for a *different* one inside every Newton iteration of every +step, which is exactly the call pattern under which allocation matters most. Measured at +`N = 128`, `p = 3`, `nq = 5`: + +| | bytes | +|:--|--:| +| `f .* q.w` temporary | 5 kB | +| `Φₐ * Diagonal(f ⊙ w)` | 47 kB | +| `(Φₐ D) * Φᵦᵀ` sparse-sparse product | 208 kB | +| **total per call** | **260 kB** | + +The `f ⊙ w` temporary is the same one `l2_projection!` no longer pays and could be removed the +same way, with the buffer the quadrature already holds. The other 255 kB are not a temporary +at all: they are the result, a freshly built `SparseMatrixCSC` with its `colptr`, `rowval` and +`nzval` allocated and its structure recomputed from scratch. + +That structure does not depend on `f`. For a fixed `(a, b)` the sparsity pattern of +`Φₐ diag(f ⊙ w) Φᵦᵀ` is the same for every coefficient — a basis function overlaps only the +`2p+1` others whose supports meet its own — so the pattern could be assembled once per +`(a, b)`, cached beside the `mixed_matrix` results, and only `nzval` refilled per +call. That turns 260 kB into zero. + +What it needs is an in-place entry point, `weighted_matrix!(A, q, f, a, b)`, since the present +signature has nowhere to write. Callers holding a matrix across steps would use it and callers +wanting a value would keep the allocating form. Deferred rather than done because it widens +the API, and because the sparse triple product would have to be written out by hand against +the cached pattern instead of delegating to `SparseArrays`, which is the part that needs to be +got right rather than merely written. From 4527bd1d1c981af94099e44f033deb244a422de9 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Fri, 21 Aug 2026 21:40:52 +0900 Subject: [PATCH 4/4] Address the review of PR #5 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `l2_projection!` took the quadrature's `f ⊙ w` buffer unconditionally, which narrowed the element types the method accepts: a complex sample threw an `InexactError` where the allocating `l2_projection` -- the very method the docstring offers as the interchangeable alternative -- returns a `ComplexF64` result, and a `BigFloat` sample was silently computed through `Float64`. The buffer is now taken only when the product lands in its element type, which is a test on types alone and so resolved at compile time: the ordinary path still allocates nothing (0 B at N = 128 on a uniform mesh), and a wider sample gets a product of its own. Verified equal to `l2_projection` to the last bit on a graded mesh. The in-place and allocating methods are separate implementations, and only the uniform mesh compared them, so the suite now checks them against each other across all three mesh families and degrees 1 to 4, and pins the complex case. The non-reentrancy warning moves from `l2_projection!` to `SplineQuadrature`, where it belongs: the shared buffer is the second piece of mutable state on that struct, not the first, the memoising `cache` behind `mixed_matrix` being equally unsynchronised. - `CompatHelper.yml` gained `julia-actions/cache` without the `actions: write` permission the action needs to prune the caches it created, which `CI.yml` grants with a comment saying so; the daily run would have logged a permission failure. - `Documenter.yml` instantiated the docs environment in a step of its own, which `julia-actions/julia-docdeploy` then did again -- it runs `Pkg.develop` and `Pkg.instantiate` on `docs/` itself. The step goes, and with it a comment describing work it was not doing. - `LICENSE.md` had lost its copyright year in the rename. - `AUTHORS.md` introduced contributors as being "in addition" to the principal developers and then listed a principal developer among them. - `.githooks/pre-push` never read the ref lines git feeds it, so it left git writing into a closed pipe and ran the whole suite for a push that only deletes a branch. Co-Authored-By: Claude Opus 5 (1M context) --- .githooks/pre-push | 22 ++++++++++++++++++++-- .github/workflows/CompatHelper.yml | 1 + .github/workflows/Documenter.yml | 9 ++------- AUTHORS.md | 3 ++- CHANGELOG.md | 11 ++++++++--- LICENSE.md | 2 +- src/quadrature.jl | 26 +++++++++++++++++++++----- test/quadrature_tests.jl | 20 ++++++++++++++++++++ 8 files changed, 75 insertions(+), 19 deletions(-) diff --git a/.githooks/pre-push b/.githooks/pre-push index b9fc049..239da8f 100755 --- a/.githooks/pre-push +++ b/.githooks/pre-push @@ -12,6 +12,24 @@ set -eu +# git feeds the hook one ` ` line per ref +# being pushed, and an all-zero local oid means that ref is being deleted. Read them: +# a push that only deletes has nothing to test, and a hook that exits without draining +# leaves git writing into a closed pipe. +refs=0 +content=0 +while read -r _local_ref local_oid _remote_ref _remote_oid; do + refs=$((refs + 1)) + case "$local_oid" in + *[!0]*) content=1 ;; + esac +done + +if [ "$refs" -gt 0 ] && [ "$content" -eq 0 ]; then + printf 'pre-push: nothing but ref deletions, skipping the test suite.\n' >&2 + exit 0 +fi + if [ "${SIMPLESPLINES_SKIP_TESTS:-0}" != "0" ]; then printf 'pre-push: SIMPLESPLINES_SKIP_TESTS is set, skipping the test suite.\n' >&2 exit 0 @@ -23,8 +41,8 @@ if ! command -v julia >/dev/null 2>&1; then exit 1 fi -# The hook runs with the working directory at the top level of the work tree already, but -# `git push` from a subdirectory of a worktree is not worth relying on that for. +# git runs its hooks from the top level of the work tree, but asking rather than assuming +# also gives the right answer in a linked worktree. root=$(git rev-parse --show-toplevel) printf 'pre-push: running the SimpleSplines test suite...\n' >&2 diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 469681a..18352eb 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -4,6 +4,7 @@ on: - cron: 0 0 * * * workflow_dispatch: permissions: + actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created contents: write pull-requests: write jobs: diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index 6902c2d..cbb6e17 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -28,13 +28,8 @@ jobs: version: '1' - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - # docs/Project.toml carries a [sources] entry for the package itself, so instantiating - # the docs environment is all that is needed to resolve it against the checkout. - - name: Instantiate the docs environment - shell: julia --project=docs --color=yes {0} - run: | - using Pkg - Pkg.instantiate() + # julia-docdeploy resolves the docs environment against the checkout and instantiates + # it itself, so there is no separate step for that here. - uses: julia-actions/julia-docdeploy@v1 env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/AUTHORS.md b/AUTHORS.md index 3489aaf..81447b3 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -13,6 +13,7 @@ Together, these two groups form "The SimpleSplines Authors" as mentioned in the ## Contributors -The following people contributed to SimpleSplines and are listed in alphabetical order: +Everyone who has contributed to SimpleSplines, the principal developers above included, in +alphabetical order: * Michael Kraus diff --git a/CHANGELOG.md b/CHANGELOG.md index 726bef8..b5a017f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -94,9 +94,14 @@ The paths a time integrator actually runs per step allocate nothing. Evaluation allocation-free at any degree and derivative order on any mesh; a `CirculantMass` solve is allocation-free; and `l2_projection!` is now allocation-free end to end on a uniform mesh, the `f ⊙ w` product going into a buffer held by the quadrature and the load vector being -formed with `mul!` straight into the output. That buffer makes `l2_projection!` -non-reentrant, which its docstring says. On a non-uniform mesh a CHOLMOD temporary remains, -CHOLMOD having no in-place `ldiv!`. +formed with `mul!` straight into the output. The buffer is taken only when the product lands +in the quadrature's element type, so a wider sample — a complex `f` — is still projected, +through a product of its own, rather than narrowed into it. On a non-uniform mesh a CHOLMOD +temporary remains, CHOLMOD having no in-place `ldiv!`. + +That buffer is mutable state on a struct that reads as immutable, as the memoising `cache` +behind `mixed_matrix` already was, so the `SplineQuadrature` docstring now warns that one +quadrature must not be shared between threads. ### Dependencies diff --git a/LICENSE.md b/LICENSE.md index 7102389..70b5dd4 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,6 +1,6 @@ MIT License -Copyright (c) The SimpleSplines Authors (see [AUTHORS.md](AUTHORS.md)) +Copyright (c) 2025-present The SimpleSplines Authors (see [AUTHORS.md](AUTHORS.md)) Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/src/quadrature.jl b/src/quadrature.jl index 465b296..184dec0 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -70,6 +70,13 @@ true support of `p+1` cells are ever nonzero and only those are computed, but keeping the array dense lets the contractions above run as one BLAS call, which is the faster arrangement at the sizes these discretisations are used at. + +!!! warning "One quadrature per thread" + A `SplineQuadrature` carries mutable state behind an otherwise read-only interface: + [`mixed_matrix`](@ref) memoises its results into a `Dict`, and [`l2_projection!`](@ref) + forms its ``f \odot w`` product in a shared buffer. Neither is synchronised, so one + quadrature must not be used from two threads at once. Give each thread its own, or stay + with the allocating [`l2_projection`](@ref), which forms its own product. """ struct SplineQuadrature{T, BT <: PeriodicBSplineBasis{T}, MO <: MassOperator{T}} basis::BT @@ -356,9 +363,12 @@ buffer held by the quadrature, the load vector is formed with `mul!` straight in the [`CirculantMass`](@ref) solve is itself allocation-free. On a non-uniform mesh the CHOLMOD solve still allocates a temporary, as [`mass_solve!`](@ref) notes. -That shared buffer is also what makes this method non-reentrant: two threads projecting -through the same [`SplineQuadrature`](@ref) at once would overwrite each other's `f ⊙ w`. -Give each thread its own quadrature, or use the allocating [`l2_projection`](@ref). +A sample whose element type is wider than the quadrature's — a complex `f` — gets its own +product instead of being narrowed into that buffer, so this method accepts exactly what +[`l2_projection`](@ref) accepts and merely stops being allocation-free there. + +The shared buffer is one of the two reasons a [`SplineQuadrature`](@ref) may not be used from +two threads at once; see the warning there. """ function l2_projection!(û::AbstractVector, q::SplineQuadrature, f::AbstractVector) length(û) == nbasis(q) || throw(DimensionMismatch( @@ -366,8 +376,14 @@ function l2_projection!(û::AbstractVector, q::SplineQuadrature, f::AbstractVect length(f) == length(q.x) || throw(DimensionMismatch( "the function was sampled at $(length(f)) points but the quadrature has " * "$(length(q.x))")) - q.scratch .= q.w .* f - mul!(û, basis_values(q, 0), q.scratch) + # The buffer has the quadrature's element type, so it can only take the product when the + # product lands in that type; a complex or extended-precision sample gets its own array + # rather than an `InexactError` or a silent narrowing. The test is on types alone, hence + # resolved when the method is compiled, so the ordinary path still allocates nothing. + S = promote_type(eltype(q.w), eltype(f)) + fw = S === eltype(q.scratch) ? q.scratch : similar(f, S) + fw .= q.w .* f + mul!(û, basis_values(q, 0), fw) mass_solve!(û, q.mass, û) return û end diff --git a/test/quadrature_tests.jl b/test/quadrature_tests.jl index 20fe773..5fe58e0 100644 --- a/test/quadrature_tests.jl +++ b/test/quadrature_tests.jl @@ -171,6 +171,26 @@ using Test l2_projection!(ŵ, q, fq) @test (@allocated l2_projection!(ŵ, q, fq)) == 0 @test ŵ ≈ û + + # The in-place and allocating methods are separate implementations -- the buffer is + # the whole point of the split -- so they are checked against each other on every + # mesh family, and not only on the uniform mesh the allocation test above needs. + for (nm, mk) in MESHES, p in 1:4 + qm = SplineQuadrature(PeriodicBSplineBasis(mk(16), p)) + fm = sin.(quadrature_nodes(qm)) + ûm = Vector{Float64}(undef, nbasis(basis(qm))) + @test l2_projection!(ûm, qm, fm) ≈ l2_projection(qm, fm) + end + + # A sample wider than the quadrature's element type gets its own f ⊙ w product + # instead of being narrowed into the shared buffer, which would be an InexactError + # here and a silent loss of precision for a BigFloat. Graded rather than uniform: + # the rfft plan of a CirculantMass takes a real argument only. + qc = SplineQuadrature(PeriodicBSplineBasis(GradedMesh(16, 2π), 3)) + fc = cis.(quadrature_nodes(qc)) + ûc = Vector{ComplexF64}(undef, nbasis(basis(qc))) + @test l2_projection!(ûc, qc, fc) ≈ l2_projection(qc, fc) + @test eltype(l2_projection(qc, fc)) == ComplexF64 end @testset "$(rpad("L2 projection converges at order p+1",76))" begin