Utilities#

backend(Wₕ::AbstractSpaceType) -> AbstractBackend

Returns the computational backend associated with the space Wₕ.

Backend{VT, MT, EP}()

Compile-time descriptor specifying vector type VT, matrix type MT, and execution policy EP.

Type parameters

• VT<:DenseVector: Concrete dense vector type (for CPU or GPU).

• MT<:AbstractMatrix: Concrete matrix type (e.g. SparseMatrixCSC{Float64, Int} or Matrix{Float64}).

• EP<:ExecutionPolicy: An ExecutionPolicy, for example Serial, Parallel, CpuPolyester or GpuKernel.

All three parameters must agree on locality: a backend is meant to be wholly host – a Vector alongside a CPU sparse or dense matrix, run under a CpuPolicy – or wholly device – once a vendor is chosen, VT, MT and the policy are all that vendor's. Any other combination (a device VT with a host MT, a device VT with a CpuPolicy, or a host VT with a GpuPolicy) is rejected at construction rather than accepted and left to fail on first use.

Throws

• ArgumentError: locality(VT) != locality(EP()).

• ArgumentError: locality(MT) != locality(VT).

• ArgumentError: EP is a GpuOffload whose wrapped device backend only supports an element type other than eltype(VT).

See also: backend, vector_type, matrix_type, execution_policy, locality.

Locality

Abstract supertype for where an array's storage lives: HostLocality for memory a CPU loop can index element by element, DeviceLocality for memory an accelerator schedules against instead.

Never declared directly – locality derives it from an array type, an ExecutionPolicy, or a Backend. So a Backend's storage and its execution policy answer the same question and can be checked against each other, rather than each independently claiming a locality that need not agree.

See also: locality, HostLocality, DeviceLocality.

HostLocality() <: Locality

Storage a CPU loop can index element by element: Vector, Matrix, SparseMatrixCSC, and, by default, any array type locality has not been taught otherwise about.

See also: Locality, DeviceLocality, locality.

DeviceLocality() <: Locality

Storage an accelerator schedules against rather than a CPU loop indexing it directly: MtlVector/MtlMatrix and any future GPU array type. Scalar indexing such storage from a host loop is the failure this trait exists to catch before it happens.

See also: Locality, HostLocality, locality.

locality(x) -> Locality

Return the Locality of x – an array type, an ExecutionPolicy, or a Backend.

Locality is derived from storage, never declared. The fallback method, locality(::Type{<:AbstractArray}) = HostLocality(), treats any array type this package has not been taught otherwise about as host memory. Any array type subtyping GPUArraysCore.AbstractGPUArray answers DeviceLocality() generically, so a GPU package extension need not add its own method for that to hold – Metal's MtlVector/MtlMatrix already subtype it, and any future CUDA/ROCm/oneAPI array type gets the right answer the same way, for free. BrambleMetalExt keeps its own MtlVector/MtlMatrix/sparse-matrix methods regardless: Julia dispatches to the more specific method, and those predate this generic one – the same idiom as ka_device: a method a package extension supplies for its own type, picked up automatically once loaded.

For an ExecutionPolicy, locality answers what the policy claims. A CpuPolicy claims HostLocality and a GpuPolicy claims DeviceLocality. For a Backend, it answers the locality of the backend's vector type VT alone, because VT is the storage a sweep actually indexes. The matrix type MT is not consulted here, since Backend's constructor checks MT against VT separately.

Examples

using Bramble
const B = Bramble
B.locality(Vector{Float64}) === B.HostLocality() &&
    B.locality(B.CpuSerial()) === B.HostLocality() &&
    B.locality(B.GpuKernel()) === B.DeviceLocality()

# output
true

See also: Locality, HostLocality, DeviceLocality, ka_device.

ExecutionPolicy

Abstract supertype for backend execution policies.

Two regimes sit under it: CpuPolicy for work a CPU loop drives, GpuPolicy for work a device schedules. A policy says where the work runs: CpuSerial walks the array element by element on one CPU thread, the one thing a device array refuses.

CpuPolicy <: ExecutionPolicy

Abstract supertype for the policies a CPU loop executes: CpuSerial, CpuThreaded and CpuPolyester.

The CPU sweeps dispatch on this rather than on each concrete policy, so a GpuPolicy reaching one is a method error with a message rather than a scalar-indexing failure several frames in.

CpuSerial() <: CpuPolicy

Sequential execution policy.

Directs grid operations and form assembly to execute via single-threaded loops. This is the default execution policy.

Also spelled Serial(): const Serial = CpuSerial.

The grid sizes at which CpuThreaded and CpuPolyester start beating this policy were measured per workload rather than assumed – see their own docstrings for the numbers.

See also: CpuThreaded, ExecutionPolicy.

CpuThreaded() <: CpuPolicy

Multithreaded execution policy.

Directs grid operations and form assembly to execute across CPU threads via static partitioning. Execution is unconditional: no per-call size thresholds are imposed. For workloads dominated by small, frequently repeated calls, use CpuSerial.

Also spelled Parallel(): const Parallel = CpuThreaded.

A call made from inside another threaded region, such as the body of a user's own Threads.@threads loop, or from a spawned task while another task's threaded region runs, executes serially on the calling task instead of throwing, and returns exactly what the same call returns at top level.

Base.Threads.@threads is the primitive, and naming it that way leaves room for the others that are not this one – Polyester's @batch and MPI. CpuPolyester is that Polyester-backed sibling: the policy type ships here, but the sweeps it selects live in the BramblePolyesterExt package extension, and requesting one without using Polyester errors the way metal_backend does without using Metal.

Where "small" ends differs widely between workloads (restriction, cell averaging, the innerₕ/normₕ reductions, D₋ₓ!, broadcasts), so a size from one does not transfer to another, and it depends on the machine. benchmark/policy_crossover.jl measures it per workload. innerₕ/normₕ (_dot/_dot_masked, src/utils/linear_algebra.jl), D₋ₓ! and the warmed broadcast axpy vₕ .= a .* uₕ .+ wₕ into a VectorElement thread under this policy.

See also: CpuSerial, CpuPolyester, ExecutionPolicy.

CpuPolyester() <: CpuPolicy

Polyester-batched execution policy.

Named for the provider, matching CpuSerial and CpuThreaded, so the name makes the using Polyester requirement self-evident.

Directs grid operations and form assembly through Polyester.jl's @batch, a primitive whose per-call overhead is low enough to pay off on grids where CpuThreaded's Threads.@threads does not. The sweeps this policy selects are implemented in the BramblePolyesterExt package extension, not here: using Polyester must be loaded before this policy reaches one of them, or the call errors naming the package, matching metal_backend's precedent. Parallel() means Threads.@threads.

Its crossover against CpuSerial falls at smaller sizes than CpuThreaded's for every workload measured, D₋ₓ! and the warmed broadcast axpy vₕ .= a .* uₕ .+ wₕ included. benchmark/policy_crossover.jl measures both on your machine.

Its reductions, innerₕ, normₕ and inner₊ with or without markers, split the sum into Threads.nthreads() fixed bands and add the band sums with sum, exactly as CpuThreaded does, and return CpuThreaded's value bitwise. The value depends only on the data and Threads.nthreads(), never on how many threads are free, so a reduction nested in another sweep, as in Rₕ!(uₕ, x -> innerₕ(vₕ, vₕ) * x[1]), returns exactly the top-level value. It agrees with CpuSerial to rounding, since the bands add the terms in another order, and on one thread it is CpuSerial's value bitwise. A task's first reduction allocates its band sums once.

A warm call under this policy allocates 0 bytes, on any grid, when the storage is a Vector or a contiguous view of one; storage that is not, such as a strided view, allocates on every call. Differences, shifts, averages, divergence, curl, εₕ!, Rₕ!, avgₕ! (masked and composite too), innerₕ, inner₊, the broadcasts (a 0-dimensional array leaf too), the weight build, linear and bilinear assembly and refill, the matrix-free product, a multigrid V-cycle, the explicit right-hand side (SemidiscretizeRHS) and the KroneckerLinearOperator product all allocate 0 B. Polyester copies what a loop captures into an argument box, and the box stays on the stack when it holds only plain arrays and isbits values. So the loops never capture a form, a space or a mesh: Bramble splits them into an isbits skeleton and their arrays, and each task rebuilds them. No call reaches Threads.@threads or Threads.@spawn.

A function you supply that captures an array, as in an avgₕ! whose source is a closure over a Vector, allocates nothing either, and sees its own objects: the loop holds a handle to a slot that stores the function, not the function. Only a function of a type beyond the 256 that have slots, a call that finds all 64 slots of its type taken (or retired by an interrupted call, for the session), and data the split cannot take apart, such as a BigFloat coefficient or a Dict, fall back to a heap-allocated box. That costs an amount per call that does not grow with the grid, and never changes the result. A task that throws rethrows on the caller with the exception type and text CpuSerial gives.

See also: CpuThreaded, CpuSerial, ExecutionPolicy.

GpuOffload{I, DB}() <: CpuPolicy
GpuOffload(device_backend::Backend, inner::CpuPolicy = CpuSerial()) -> GpuOffload

A CpuPolicy that wraps a host inner policy alongside a device Backend device_backend, for host-side Backends that offload individual calls (Rₕ!, avgₕ!) to the device rather than running them on the host.

Carries no runtime fields: I and DB are type parameters only, so a Backend{VT, MT, GpuOffload{I, DB}} still rebuilds EP() the way every other ExecutionPolicy does. locality answers HostLocality for it, as for any CpuPolicy: the Backend it configures is a host backend, and device_backend is only where the offloaded calls are dispatched, not this backend's own storage. execution_policy on a Backend configured with a GpuOffload policy returns I(), the wrapped inner policy instance, not the GpuOffload itself – every existing concrete-policy dispatch site sees the inner policy and behaves exactly as it does without a GpuOffload in play.

Throws

• ArgumentError: GpuOffload(device_backend, ...) with a device_backend whose locality is HostLocality – nothing would be offloaded.

• ArgumentError: constructing a Backend with a GpuOffload{I, DB} policy whose requested vector element type cannot be represented by device_backend's own vector element type (e.g. Float64 requested against a Metal device_backend, which only supports Float32 and Float16).

See also: CpuPolicy, Backend, metal_backend, execution_policy.

GpuPolicy <: ExecutionPolicy

Abstract supertype for the policies a device executes, currently GpuKernel.

See also: CpuPolicy, ExecutionPolicy.

GpuKernel() <: GpuPolicy

Device execution policy: work is dispatched to the accelerator rather than to a host loop.

What metal_backend carries by default. A GPU is massively parallel and cannot execute serially, so Serial() would be a false statement about the hardware and would send CPU assembly loops at device arrays.

The name states the launch mechanism only: work is dispatched to a device kernel. A kernel launch only enqueues onto the device's queue; a host read of the result, or ka_synchronize, is the barrier. There is no handle to wait on. Kernels on one queue run in its order, so a chain of operators needs no synchronise between its steps. The memory locality lives in the Locality trait, which this policy claims through locality.

See also: GpuPolicy, ExecutionPolicy.

Serial

Alias for CpuSerial, the spelling the package's call sites, tests and benchmark baseline keys use.

Parallel

Alias for CpuThreaded.

execution_policy(backend::Backend{<:Any, <:Any, EP}) -> EP
execution_policy(::Type{<:Backend{<:Any, <:Any, EP}}) -> EP

Return the ExecutionPolicy instance configured for backend.

For a GpuOffload backend this is the wrapped inner policy I(), not the GpuOffload wrapper, so every concrete-policy dispatch site sees the policy that runs the host work.

execution_policy(Ωₕ::AbstractMeshType) -> ExecutionPolicy

Return the ExecutionPolicy (Serial or Parallel) of the Backend associated with mesh Ωₕ.

vector(backend::Backend{VT}, n::Integer) -> VT

Allocate an uninitialized vector of length n using vector type VT configured in backend.

Arguments

• backend: Target backend instance.

• n: Number of vector elements.

Which constructor is used is decided by supports_undef_construction(VT), from the type alone.

Throws

• ErrorException: If VT does not define the constructor its trait declares.

matrix(backend::Backend{<:Any, MT}, n::Integer, m::Integer) -> MT

Allocate a matrix of dimensions n × m using matrix type MT configured in backend.

For dense matrix types, allocates uninitialized storage via MT(undef, n, m). For sparse matrix types (SparseMatrixCSC), allocates an empty sparse matrix via spzeros(T, Ti, n, m).

Arguments

• backend: Target backend instance.

• n: Number of rows.

• m: Number of columns.

Which constructor is used is decided by supports_undef_construction(MT), from the type alone.

Throws

• ErrorException: If MT does not define the constructor its trait declares.

vector_type(backend::Backend{VT}) -> Type{VT}
vector_type(::Type{<:Backend{VT}}) -> Type{VT}

Return the vector type VT configured for backend.

matrix_type(backend::Backend{<:Any, MT}) -> Type{MT}
matrix_type(::Type{<:Backend{<:Any, MT}}) -> Type{MT}

Return the matrix type MT configured for backend.

backend_types(backend::Backend{VT, MT, EP}) -> Tuple{Type, Type{VT}, Type{MT}, Type{Backend{VT, MT, EP}}}
backend_types(::Type{<:Backend{VT, MT, EP}}) -> Tuple{Type, Type{VT}, Type{MT}, Type{Backend{VT, MT, EP}}}

Return a 4-tuple containing (eltype(VT), VT, MT, Backend{VT, MT, EP}).

supports_undef_construction(::Type{AT}) -> Bool

Whether an array type builds from an undef initializer, AT(undef, dims...), rather than from its dimensions alone, AT(dims...).

The contract a custom array type opts into to be usable as a backend's vector_type or matrix_type. true by default, which is what every AbstractArray in Base and in every GPU package this repository has seen answers; a type constructed from its size alone declares otherwise:

Bramble.supports_undef_construction(::Type{<:MySizedArray}) = false

Read at compile time from the type, so the branch it guards folds away and the compiler sees through the allocation path. The three shipped backends never reach it – Vector, DenseMatrix and SparseMatrixCSC each have their own @inline method below – so it is the stated contract for the types that do.

backend_eye(backend::Backend, n::Integer) -> AbstractMatrix

Construct an \(n \times n\) identity matrix matching the matrix type configured in backend.

backend_zeros(backend::Backend, n::Integer) -> AbstractMatrix

Construct an \(n \times n\) zero matrix matching the matrix type configured in backend.

ka_device(be) -> KernelAbstractions.Backend

Return the KernelAbstractions.jl device backend that runs kernels for the Bramble Backend be.

This is the extension contract every device backend implements: BrambleMetalExt adds a method for a Backend built from MtlVector/MtlMatrix, returning Metal.MetalBackend(), and a future GPU extension does the same for its own array types. Every @kernel declared in BrambleKernelAbstractionsExt is written against KernelAbstractions.Backend alone, so it becomes available on a new device the moment that device supplies this one method.

Throws

• ErrorException: if no loaded package defines ka_device for be. Requires using KernelAbstractions and the package that supplies be's device (for example using Metal for a Metal backend).

ka_synchronize(x) -> Nothing

Block the calling thread until every KernelAbstractions.jl kernel and device transfer already queued against x's device backend has completed.

Under GpuKernel – the only GpuPolicy there is – a kernel launch only enqueues onto the device's own command queue and returns, so a chain of operators (D₋ₓ into D₋ᵧ, say) pipelines instead of paying a host round-trip after each step. Kernels enqueued on the same queue still run in that queue's order, so this is not needed between chained device calls, only at a genuine host boundary:

• converting a device array to a host one (Array(...), host_points, ...), or any other host-side read, reduction or assertion – though a GPUArrays reduction that returns a plain host scalar (sum, which _dot's device method in src/utils/linear_algebra.jl uses, and so innerₕ/normₕ too) already synchronises by fetching that scalar, with no call to this function needed;

• a write that reaches device memory some way other than a @kernel launch on that same queue – a plain copyto!, which queues asynchronously exactly like a kernel launch does, but is not one of the launches above. _flush_device_scatter! (src/assembly/bilinear_traversal.jl) is this function's first caller for exactly that reason: it ends a device-resident matrix's assembly with copyto!(A.nzVal, mirror.nzval), and calling this right after is what keeps assemble/assemble! from returning before that write lands (a race found at n = 513 over repeated assemblies, invisible at a small test size);

• an explicit call from user code that needs a hard barrier before doing something this package cannot see, such as timing a device computation in isolation.

Same idiom as ka_device: a helpful error naming the packages to load, not a bare MethodError, and the real method comes from BrambleKernelAbstractionsExt (synchronize(get_backend(x))) – written against KernelAbstractions alone, so a future GPU backend inherits it for free.

Throws

• ErrorException: if no loaded package defines ka_synchronize for x. Requires using KernelAbstractions.

metal_sparse_csr(A::SparseMatrixCSC) -> MetalSparseMatrixCSR

Convert a host SparseMatrixCSC to a sparse matrix in compressed sparse row (CSR) format, stored in Metal device memory.

Requires using Metal in the caller environment. The conversion is sparse-to-sparse throughout – A is never densified.

Arguments

• A: The host sparse matrix to convert.

Throws

• ErrorException: If Metal.jl is not loaded.

metal_sparse_csc(A::SparseMatrixCSC) -> MetalSparseMatrixCSC

Convert a host SparseMatrixCSC to a sparse matrix in compressed sparse column (CSC) format, stored in Metal device memory.

Requires using Metal in the caller environment. SparseMatrixCSC's own storage already is CSC, so no format conversion happens – its fields are placed in device memory as-is.

Arguments

• A: The host sparse matrix to convert.

Throws

• ErrorException: If Metal.jl is not loaded.

gpu_backend(::Type{T} = Float32; policy::ExecutionPolicy = GpuKernel()) -> Backend

Construct a GPU Backend for whichever accelerator extension is loaded and has a functional device.

Detects a loaded GPU package extension with Base.get_extension, checks _gpu_functional for that device, and forwards to the backend's own constructor: currently metal_backend, under using Metal with Metal.functional() true. A CUDA or AMDGPU extension joins this dispatch once it exists.

Arguments

• T: Floating-point element type (default: Float32).

Keywords

• policy: Execution policy instance (default: GpuKernel).

Throws

• ErrorException: two distinct diagnoses, deliberately not sharing a message. If no GPU extension is loaded at all, the message names the package to load, chosen from the host architecture: Metal on Apple Silicon, CUDA on Linux/Windows, and a generic "no GPU hardware" message otherwise. If a GPU extension is loaded but its device is not functional (e.g. using Metal on a host without working Metal drivers), the message says so instead – that is a driver, hardware or virtualisation problem, not a missing using.

See also: metal_backend, backend.

metal_backend(::Type{T} = Float32; policy::ExecutionPolicy = GpuKernel()) -> Backend

Construct a Metal GPU Backend backed by Metal.jl arrays.

Requires using Metal in the caller environment. Apple Silicon GPUs support Float32 and Float16, but do not support 64-bit floating point arithmetic.

A function projected on this backend has to be GPU-compilable. Rₕ(W, f) on a metal_backend space compiles f into a device kernel rather than calling it from a host loop, so f must be device-compilable: no Float64 literals (a bare literal like 0.5 forces double precision, which Apple Silicon GPUs do not support – write x[1] / 2, not x[1] * 0.5), no allocations, and no calls to a non-inlineable or host-only function. x -> sin(x[1]) compiles; anything capturing a boxed value or calling out to the host does not, and fails at compile time with a dynamic-invocation error rather than at the call site.

Arguments

• T: Floating-point element type (Float32 or Float16, default: Float32).

Keywords

• policy: Execution policy instance (default: GpuKernel). A CpuPolicy is rejected at construction, not accepted and left to fail on first use: Metal's device arrays answer DeviceLocality to locality, a CpuPolicy answers HostLocality, and Backend's inner constructor requires the two to agree.

Throws

• ErrorException: If Metal.jl is not loaded.

• ArgumentError: If policy is a CpuPolicy (locality mismatch), or if T is a numeric type other than Float16 and Float32 (with Metal.jl loaded), such as Float64.

csr_backend(::Type{T} = Float64; policy::ExecutionPolicy = Serial()) -> Backend

Construct a SparseMatrixCSR Backend backed by SparseMatricesCSR.jl.

Requires using SparseMatricesCSR in the caller environment. A finite-difference stencil is assembled row by row, which compressed sparse row storage reaches without the column scatter a SparseMatrixCSC assembly needs.

Arguments

• T: Coordinate and scalar element type (default: Float64).

Keywords

• policy: Execution policy instance (default: Serial).

Throws

• ErrorException: If SparseMatricesCSR.jl is not loaded.

profile_backends() -> BackendProfile

Time one sweep kernel under each execution policy this session can run on the host, for sizes 2^10, 2^12, ..., 2^22, and report where each policy overtakes Serial().

Each policy sweeps a vector through the same loop the package's own weight and scatter builds use, with an isbits kernel (i -> sqrt(i) * a + b): one warm run, then the best of three runs. The candidates are Serial(), Parallel(), and CpuPolyester() once using Polyester has loaded its extension. With using Metal on a functional device a GpuKernel() row is added, timed with a device synchronisation. Every row is timed on the same element type, so the rows compare: Float64, or Float32 for all of them once the Metal row is present (the only element type Metal offers), in which case the figures are not drop-in figures for a Float64 problem. The header of the table names the type. The host policies take well under a second including compilation. With Metal loaded, the first call is slow because it is often the session's first GPU operation, which compiles Metal's own code as well as the sweep kernel; loading Polyester before Metal makes that first GPU operation markedly slower (an upstream interaction), and loading it after Metal reduces it. Later calls take well under a second. It is meant to be called by hand, once per session, and is not run by backend or gridspace.

The result prints as a table of times, the crossover of each policy against Serial() (the first size from which it stays at least 1.2 times faster to the end of the sweep, or never) and the backend(policy = ...) spelling to use. The table times one sweep kernel, not assembly: a policy that wins here is a good first choice, not a guarantee for your form. The times are plain data in the eltype, sizes, labels and times fields.

Returns

• BackendProfile: eltype, sizes, policy labels, spellings and the times matrix in seconds.

Examples

using Bramble
Bramble.profile_backends()

Matrix storage

profile_backends times execution policies only. It does not compare csr_backend with the default SparseMatrixCSC storage: the matrix type only matters during assembly and solves, and how much it matters depends on the form, the boundary conditions and the solver. To compare the two on your own problem, load SparseMatricesCSR, build the same mesh, space and form under each backend, and time assemble, assemble! and the solve you use:

using Bramble, SparseMatricesCSR, BenchmarkTools

Ω = domain(box((0.0, 0.0), (1.0, 1.0)))
for be in (backend(), csr_backend())
    Wₕ = gridspace(mesh(Ω, (257, 257); backend = be))
    a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))
    A = assemble(a; dirichlet = :boundary)
    @btime assemble!($A, $a; dirichlet = :boundary)
end

Use the policy profile_backends recommended in both backends, e.g. csr_backend(; policy = Parallel()), so the comparison isolates the storage.

See also: backend, csr_backend, Parallel.

BackendProfile

Result of profile_backends: plain data, printed as a table.

Fields

• eltype::DataType: Element type every policy was timed on.

• sizes::Vector{Int}: Number of points of each sweep.

• labels::Vector{String}: One label per policy; the first is the Serial() baseline.

• spellings::Vector{String}: The expression that selects each policy.

• times::Matrix{Float64}: Seconds, times[i, j] for sizes[i] under policy j.

sparse_refactor!(fact::Factorization, A::SparseMatrixCSC) -> Factorization
sparse_refactor!(fact::Factorization, a::BilinearForm; kwargs...) -> Factorization

Alias for refactor!. Recomputes numerical values of fact reusing the existing symbolic factorization.

Bramble.PRECOMPILE_WORKLOAD -> Bool

Whether the @compile_workload blocks – this file's own, and each package extension's (ext/*.jl, gpena/Bramble.jl#196) – run at precompile time. Read from the "precompile_workload" Preferences.jl key, defaulting to true.

using Preferences, Bramble
set_preferences!(Bramble, "precompile_workload" => false)

disables every workload, core and extensions alike, since they all gate on this same constant; the change takes effect on the next load, no manual cache clearing needed.

public, not exported: a real, cross-module contract every extension's own precompile gate reads, but not something an ordinary caller uses.