Solvers and preconditioners#

pde_solve(A::SparseMatrixCSC, F::AbstractVector; solver = :default, sym = :auto, mean_zero = nothing, gauge = nothing, kwargs...) -> Vector
pde_solve(fact::Factorization, F::AbstractVector) -> Vector

Solve A u = F (or fact u = F) and return u.

With no keyword arguments this is A \ F, unless AppleAccelerate.jl is loaded on macOS and the system is symmetric (see :default below), in which case it is accelerate_solve(A, F; sym). The name exists so that reverse-mode AD tools have one function to attach an adjoint rule to, which is what makes a whole θ -> assemble -> pde_solve -> J(u) chain differentiable.

Keywords

• solver: :default – A \ F on Linux, Windows, macOS without AppleAccelerate.jl loaded, and macOS on an unsymmetric A; on macOS with using AppleAccelerate in effect and A symmetric, dispatches to accelerate_solve instead (same as passing solver = :accelerate explicitly). Symmetry is issymmetric(A) under sym = :auto (the default), or trusted outright from an explicit sym = :spd/:definite/:symmetric (and conversely :unsymmetric skips straight to A \ F) – the caller has already asserted the property, so it is not checked again. The narrowing exists because Accelerate is only a win on the symmetric factorisations it reaches (SPD/Cholesky, LDLᵀ). On forms assembled on 2D grid spaces it won on a symmetric Poisson-plus-mass system and lost on an unsymmetric convection-diffusion one. This only ever narrows which solver runs on macOS; Linux and Windows are never affected, and a matrix Accelerate cannot factor (e.g. non-square) throws from accelerate_solve exactly as solver = :accelerate would – :default never silently falls back to \ after picking Accelerate. It picks Accelerate only for Float32, Float64 and their complex types (the element type of A and of the solution), so an integer A goes to A \ F, or in the mean_zero/gauge solve to lu, which converts it to Float64, and a type that backend cannot factor (BigFloat) throws an ArgumentError naming it. :suitesparse (CHOLMOD/UMFPACK), :spqr (sparse QR, for a least-squares or rectangular A), :accelerate (Apple libSparse, needs AppleAccelerate.jl, honoured unconditionally regardless of symmetry), :mumps (needs MUMPS.jl), :sparspak (pure Julia, needs Sparspak.jl).

• kwargs: passed to the chosen backend (for :accelerate and :mumps the keywords of accelerate_solve and mumps_solve), which throws on one it does not take. :default takes only kind and :sparspak takes none. Both throw an ArgumentError on any other keyword. kind and sym are read only when the Accelerate route runs. On A \ F they are accepted and ignored, so kind = :bogus throws only where Accelerate is taken.

• sym: symmetry hint for :suitesparse, :accelerate and :mumps (and, on macOS, for :default's own choice of solver – see above). :auto (default) detects it; :spd/:definite/1, :symmetric/2 and :unsymmetric/0 state it.

• mean_zero: a weight vector w (for instance the assembled innerₕ(1, v)), or nothing (default). With w, pde_solve solves the bordered system [A w; wᵀ 0] [u; λ] = [F; 0] and returns u, the solution with discrete zero average \(w^T u = 0\) of a pure-Neumann system whose A annihilates constants, of the promoted float type of A, F and w; A is not modified. The bordered matrix is indefinite, so sym is not used. For a real symmetric A and real w, :default factors it by LU (or by Accelerate's pivoting \(LDL^T\)) and :accelerate and :mumps by \(LDL^T\); for any other A they use LU, as do :suitesparse (UMFPACK) and :sparspak. :spqr throws an ArgumentError, and a w, F or A of mismatched size a DimensionMismatch. Without mean_zero or gauge, a square A with ‖A 1‖∞ ≤ 1e-10 ‖A‖∞ (constants in its kernel) logs a warning naming mean_zero, and the solve proceeds.

• gauge: a matrix Z with one column per kernel vector to fix (for instance Bramble.ghost_kernel(A, Wₕ), the ghost-only kernel of a form such as innerₕ(Δₕ(u), Δₕ(v)) on a ghost space), or nothing (default). With Z, pde_solve solves [A Z; Zᵀ 0] [u; μ] = [F; 0] and returns u, the solution with \(Z^T u = 0\), which makes its component in the kernel deterministic. With mean_zero = w too, the border is [w Z]. The routing, the element type and the refusals are those of mean_zero, a Z whose row count is not that of A being a DimensionMismatch, except that :sparspak, which factors without pivoting, throws an ArgumentError. A complex A or Z is accepted; the border is transposed, not conjugated. Z is used as given, and the values in the kernel are deterministic only when its columns span the whole kernel. A Z that spans only part of it is not detected: the factorisation usually proceeds (a SingularException is a matter of rounding), the physical part is still the solution of the quotient system, and the kernel components Z leaves out come back arbitrary and different from one backend to another. The solution is only as accurate as the conditioning of A away from its kernel allows: on the forms that produce a ghost-only kernel, near 1e-8 relative.

Returns

• Vector: the solution, of the promoted element type of A and F.

A SparseMatrixCSR (SparseMatricesCSR.jl) A is also accepted: converted to SparseMatrixCSC first (see docs/src/internals/csr_solvers.md), then solved exactly as above.

Reverse-mode differentiation

With ChainRulesCore.jl loaded, BrambleChainRulesExt's rrule solves the adjoint system Aᴴ λ = ∂J/∂u once and returns ∂J/∂A = -λ u' (conjugate transpose, the ChainRules convention for a complex A; -λ uᵀ for real data) restricted to A's sparsity (never densified) and ∂J/∂F = λ. Without keywords the adjoint solve reuses the forward solve's own factorisation. The rule accepts the keywords of pde_solve: it forwards them to the primal solve and to the adjoint solve of sparse(A'), so solver and sym apply to Aᴴ as well, and a sym hint must hold for Aᴴ too. A SparseMatrixCSR A has no rule (ChainRules sees a method with no rrule, and differentiates the conversion and the solver's internals). using Enzyme is enough for the same adjoint through BrambleEnzymeExt's native EnzymeRules rule. It takes a SparseMatrixCSC A, a vector F and the keywords of pde_solve, forwarded as the rrule forwards them, and differentiates in reverse mode at width 1: a batched reverse throws an ArgumentError. It conjugates u as the rrule does. Do not call Enzyme.@import_rrule, whose bridge drops the cotangent's explicit zeros from nzval and returns a wrong gradient without any error.

A closure capturing a grid space or a form needs Enzyme.Const(f) and Enzyme.set_runtime_activity(Enzyme.Reverse), as the automatic differentiation tutorial describes for every other path. A gradient with respect to a runtime scalar coefficient of the form should be a Float64 or a Ref: an Integer coefficient known only at run time costs form its inferred return type, which Enzyme's type analysis rejects.

Forward-mode differentiation

With ForwardDiff.jl loaded, a Dual matrix or load of the plain solve is solved by BrambleForwardDiffExt, which factors the Float64 values of A once, by the route the Float64 solve of the same keywords takes, and solves every partial against that factor: with \(A₀ = value.(A)\) and \(b₀ = value.(F)\), the value is \(x₀ = A₀⁻¹b₀\) and the partial \(k\) is \(xₖ = A₀⁻¹(bₖ − Aₖx₀)\). A nested Dual recurses through the same factor. A non-square Dual matrix throws an ArgumentError: a least-squares solve has a residual term that this derivative omits.

The bordered solve (mean_zero, gauge, and the rank-one term of a Bramble.RankOneUpdate) differentiates the same way on the bordered matrix: the primal bordered system of the values is factored once, by the route the Float64 solve takes, and each partial solves its bordered right-hand side \([Fₖ − Aₖu₀; 0]\) against that factor. w and gauge's Z are constants: a Dual one throws an ArgumentError naming it, as does a RankOneUpdate over a Dual matrix, whose σ takes the matrix's element type. pde_solve(fact, F) is unchanged, so with a Dual load it takes the route it always took.

Examples

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
u = pde_solve(A, F)

using SuiteSparse
u_spd = pde_solve(A, F; solver = :suitesparse, sym = :spd)

# pure Neumann: fix the constant by a discrete zero average
w = assemble(form(W, v -> innerₕ(x -> 1.0, v)))
u0 = pde_solve(A_neumann, F_neumann; mean_zero = w)

# a ghost-only kernel: fix the ghost values no row of A sees
Wg = gridspace(Ωₕ; ghosts = :boundary)
Ag = assemble(form(Wg, Wg, (u, v) -> innerₕ(Δₕ(u), Δₕ(v)) + innerₕ(u, v)))
ug = pde_solve(Ag, Fg; gauge = Bramble.ghost_kernel(Ag, Wg))

For the Bramble.RankOneUpdate that a form with a rank-one term assembles to, see pde_solve(K::Bramble.RankOneUpdate, F) below.

See also assemble, sparse_factorize, suitesparse_solve, suitesparse_qr_solve, accelerate_solve, mumps_solve, sparspak_solve, linear_problem.

pde_solve(K::Bramble.RankOneUpdate, F::AbstractVector; solver = :default, sym = :auto, mean_zero = nothing, gauge = nothing, kwargs...) -> Vector

Solve (A + σ wₛ wₜᵀ) u = F for the Bramble.RankOneUpdate K that assemble returns for a form holding a rank-one term σ * innerₕ(g, u) * innerₕ(h, v), and return u. The dense matrix is never formed: pde_solve solves the sparse bordered system

whose second row gives μ = σ wₜᵀu, so its first row is A u + σ wₛ (wₜᵀu) = F. This holds for any A, so a pure-Neumann A, singular alone, is solved when the penalty fixes its constant. For data with 1ᵀF = 0 the solution of A + σ w wᵀ (wₛ = wₜ = w, with 1ᵀA = 0) has wᵀu = 0, so in exact arithmetic it is pde_solve(A, F; mean_zero = w) for every σ ≠ 0. In floating point a small σ costs accuracy, since −1/σ then dwarfs A's entries: with σ = γ / hₘₐₓ the relative difference from the mean_zero solve is about 1e-14 / γ on the meshes tested (1e-15 for γ ≥ 1, 1e-8 at γ = 1e-6, 1e-2 at γ = 1e-12). Choose γ ≥ 1e-2 for ten digits; a large σ costs nothing. K, K', transpose(K) and transpose(K') are all accepted.

σ is K.σ, the value read when K was last assembled: after changing a Ref factor of σ, call assemble!(K, a) before solving.

Keywords

The keywords and backends of pde_solve(A, F; mean_zero), which this solve shares, its routing and element type included: solver = :default factors by LU, or by Accelerate's \(LDL^T\) when the bordered matrix is symmetric (A real symmetric and wₛ == wₜ); :suitesparse and :sparspak by LU, :accelerate and :mumps by \(LDL^T\) or LU. sym describes A and is not used. :spqr throws an ArgumentError. mean_zero = w and gauge = Z add their columns to the border after wₛ, as they do for a matrix.

Edge cases

• With σ = 0 the term vanishes, and pde_solve(K.A, F; ...) runs with every keyword, so the result is that of the plain matrix (including its warning for a pure-Neumann A).

• A negative or complex σ solves the same bordered system, valid whenever it is nonsingular. A complex σ or vector makes the border complex, transposed and not conjugated. For a nonsingular A a negative σ with 1 + σ wₜᵀA⁻¹wₛ = 0 makes A + σ wₛ wₜᵀ singular; the factorisation does not notice, and the solution comes back huge without an error. A pure-Neumann A with 1ᵀwₛ and 1ᵀwₜ nonzero is nonsingular for every σ ≠ 0.

• With wₛ ≠ wₜ, a non-symmetric term, the border [A wₛ; wₜᵀ −1/σ] is routed as unsymmetric.

• A Dirichlet-modified K from assemble(a; dirichlet), with wₛ zero at the Dirichlet rows, keeps those rows as identity rows of the bordered matrix; symmetrize!'s K solves as well.

• A singular A the term cannot fix: when ‖A 1‖∞ ≤ 8 ε ‖A‖∞ (constants in its kernel, as for pure Neumann conditions) and |1ᵀwₜ| ≤ 8 ε ‖wₜ‖₁ (the trial functional does not see constants), or the same for 1ᵀA and wₛ, A + σ wₛ wₜᵀ is singular and an ArgumentError names the fix (a functional with nonzero sum, or mean_zero). ε is the machine epsilon of A's element type. The check runs only without mean_zero and gauge; any other singular case reaches the factorisation.

• A K.A that is a SparseMatrixCSR, or the adjoint or transpose of a sparse matrix, is converted to SparseMatrixCSC; any other non-sparse K.A throws an ArgumentError.

Returns

• Vector: the solution, of the float type promoted from K and F (for σ = 0, the plain matrix's result type). F is one vector, as for a matrix.

Examples

a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)) + σ * innerₕ(1, u) * innerₕ(1, v))
K = assemble(a)
u = pde_solve(K, F)  # pde_solve(A, F; mean_zero = w) when sum(F) == 0

ghost_kernel(A, Wₕ::GhostGridSpace; rtol = 64 eps) -> SparseMatrixCSC

A sparse ndofs(Wₕ) × k matrix whose columns span the ghost-only kernel of A, a matrix assembled on the ghost space Wₕ (built by gridspace(Ωₕ; ghosts = ...) or constrained): the vectors z with A z = 0 that are zero at every physical degree of freedom.

The kernel is computed one physical node at a time. The ghosts at a node that remain unknowns, one for each ghosted facet through it without a relation, give m columns of A; restricted to their nonzero rows they form a small dense block, and each of its right singular vectors whose singular value is at most rtol times the larger of the block's largest one and the largest column 2-norm of A becomes a column of the result, with at most m entries. The second scale makes a ghost column that is zero only to rounding (a conormal K that cancels a coefficient) a kernel vector, as the dense spectrum calls it. The columns have unit 2-norm and their first nonzero entry is real and positive; the nodes come in the mesh's linear order, so the result depends on A and Wₕ only.

For a form that reads the ghosts at a corner through Δₕ alone, such as innerₕ(Δₕ(u), Δₕ(v)) + innerₕ(u, v), the corner's ghost columns are parallel, so each 2D corner gives one column and each 3D corner two, and each node inside a 3D edge one: k = 4 in 2D, k = 4 Σₑ (nₑ - 2) + 16 in 3D, with nₑ the points along axis e. With a single ghosted facet no node carries two ghosts, and k = 0 for that form. A form that never reads a ghost leaves it in the kernel, one column each.

A kernel vector spread over the ghosts of several nodes is not found. Such vectors arise when the form differences the ghosts along their own facet: with ghosts = (:xmin,), innerₕ(∇ₕ(Δₕ(u))[2], v) annihilates the vector that is constant on the :xmin ghosts, and ghost_kernel returns no column for it. A gauged solve with such a Z is not refused and usually factors: the physical part is still right, and the components left out are arbitrary.

ghost_kernel is public, not exported: write Bramble.ghost_kernel.

Throws an ArgumentError when A is not an AbstractMatrix of size ndofs(Wₕ) × ndofs(Wₕ), or when Wₕ is not a GhostGridSpace. A space whose relations eliminate every ghost has no ghost-only kernel, and the result has no columns.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 2.0)), (6, 7), (false, false))
Wₕ = gridspace(Ωₕ; ghosts = :boundary)
A = assemble(form(Wₕ, Wₕ, (u, v) -> innerₕ(Δₕ(u), Δₕ(v)) + innerₕ(u, v)))
Z = Bramble.ghost_kernel(A, Wₕ)    # 68 × 4, two entries per column, one column per corner
norm(A * Z)                        # rounding

neumann_compatibility(problem, Ωₕ[, domain_markers]) -> NamedTuple

Check that the data of a pure-Neumann problem are compatible, by building the problem on the mesh Ωₕ and on copies of it refined once and twice at the midpoints. problem(Ωₕ) returns (A, F, w): the assembled matrix A, whose kernel holds the constants, the load F and the weights w of the discrete average (for instance the assembled innerₕ(1, v)). Ωₕ itself is not changed. The copies are refined by iterative_refinement!: pass the domain's domain_markers (markers(Ω)) when Ωₕ carries markers of its own, which the refinement re-evaluates on the new points; without them such a mesh is refused.

On each mesh the system is solved with pde_solve(A, F; mean_zero = w), and the bordered solve's multiplier \(\lambda = w^\top (F - A u) / w^\top w\) measures the part of F that no u reaches. The signed load it removes and its relative defect are

both independent of the scale of w (\(\delta\) is zero for zero data). For compatible data \(\mu\) is a quadrature error and falls under refinement, by about 4 when the spacings halve; an incompatible datum, such as a source with a nonzero mean added to compatible flux data, leaves a \(\mu\) that tends to the datum's integral and keeps its sign.

The data are reported incompatible when the defect \(\delta_f\) on the refined mesh is above 5e-3 and the ratio \(\mu_f / \mu_{ff}\) of the removed loads on the once- and twice-refined meshes lies in \([0.4, 2.5]\); Ωₕ itself does not decide, since its quadrature error can be pre-asymptotic. The thresholds come from measurement on random, graded and one-axis-uniform 2D meshes of the unit square and of \([0, 1] \times [0, 2]\), with 9 to 129 points per axis: the ratio was negative or at least 3.39 for compatible harmonic data, 0.99 to 1.01 for a source shifted by 0.2, and 0.64 to 1.9 for a shift of 5 or ±20 against the larger quadrature error of \(e^{5x} \cos 5y\) fluxes (random meshes of the square, 9 to 33 points).

The check has known limits:

• An incompatibility whose refined defect is below 5e-3, or which is small against the quadrature error of the refined meshes (seen on coarse meshes of the rectangle), is not detected.

• Compatible data whose flux is singular at a boundary point (as for \(u = \operatorname{Re} z^{\alpha}\) near a corner, \(\alpha < 1\)) can be reported incompatible: their quadrature error falls like \(h^{\alpha}\), a ratio of \(2^{\alpha} < 2\) per refinement that no ratio of two or three meshes tells from an incompatibility.

• When \(\delta_f\) is above 5e-3 the problem is built and solved on the twice-refined mesh, with 16 times the points of Ωₕ: at 129 points per axis that took 20.8 s and 3.4 GiB, against 1.4 s when the floor decides.

An incompatible datum logs a warning reporting the relative defects, the multipliers and the ratios; nothing is thrown. The result is the NamedTuple (; defect, refined_defect, λ, refined_λ, ratio, second_ratio, compatible), with ratio = μ / μ_f, the removed loads on Ωₕ and its refinement (infinite when μ_f is zero), second_ratio = μ_f / μ_ff, which decides, or nothing when the floor decided and the second refinement was not built, and compatible the negation of the rule above.

neumann_compatibility is public, not exported: write Bramble.neumann_compatibility.

Throws an ArgumentError when A, on any of the meshes, does not annihilate the constants (\(\| A \mathbf{1} \|_\infty > 1000\, \varepsilon \| A \|_\infty\), with \(\varepsilon\) the machine epsilon of A's element type), and the errors of iterative_refinement! for a mesh it cannot refine.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (17, 17), (false, false))
function problem(Ωₕ)   # -Δu = f + s with ∂ₙu = ∂ₙ(cos πx cos πy): compatible only for s = 0
    Wₕ = constrained(gridspace(Ωₕ; ghosts = :boundary),
        u -> (:xmin => ∂ₙ(u) => 0.0, :xmax => ∂ₙ(u) => 0.0,
              :ymin => ∂ₙ(u) => 0.0, :ymax => ∂ₙ(u) => 0.0))
    f = x -> 2π^2 * cos(π * x[1]) * cos(π * x[2]) + 0.2
    A, F = assemble(form(Wₕ, Wₕ, (u, v) -> innerₕ(-Δₕ(u), v)), form(Wₕ, v -> innerₕ(f, v)))
    w = assemble(form(gridspace(Ωₕ), v -> innerₕ(x -> 1.0, v)))
    return A, F, w
end
r = Bramble.neumann_compatibility(problem, Ωₕ)   # warns: λ ≈ 0.2 on every mesh
r.compatible                                    # false

amg_preconditioner(A::AbstractMatrix; method = :smoothed_aggregation, kwargs...) -> MultiLevel
amg_preconditioner(a::Union{BilinearForm, Bramble.GhostBilinearForm}; method = :smoothed_aggregation,
                    dirichlet = nothing, dirichlet_components = nothing, kwargs...) -> MultiLevel

Build an algebraic multigrid hierarchy for the symmetric positive-definite matrix A – or for a assembled with the given Dirichlet conditions – ready to turn into a preconditioner with AlgebraicMultigrid.aspreconditioner.

The elliptic SBP forms Bramble assembles (Poisson, variable-coefficient diffusion, Helmholtz) have condition numbers scaling as O(h^-2): a direct sparse solve suffers severe fill-in past a few hundred thousand degrees of freedom, and an unpreconditioned Krylov method needs O(h^-1) iterations. Algebraic multigrid builds its coarse-grid hierarchy from the graph of A alone, giving a preconditioned Krylov method grid-independent, O(1) iteration counts instead.

Keywords

• method: :smoothed_aggregation (default) or :ruge_stuben, AlgebraicMultigrid's two hierarchy constructions. The two return different MultiLevel types, so a method known only at run time makes the result type a Union. A literal method, or the default, gives a concrete type.

• Every other keyword forwards to the chosen AlgebraicMultigrid constructor.

The BilinearForm method matches assemble(a::BilinearForm; ...)'s own default and does not symmetrize – Dirichlet rows become eₖ, but the matching columns are left alone, so A is not exactly symmetric even though the underlying operator is. AMG still builds a usable hierarchy from it, but for the SPD matrix the theory assumes, assemble through assemble(a, l; symmetrize = true) and pass that matrix to the AbstractMatrix method instead.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

A form on ghost spaces (GhostGridSpace) is assembled to its reduced matrix, assemble(a; dirichlet, dirichlet_components), and coarsened as any other.

Requires AlgebraicMultigrid.jl; call using AlgebraicMultigrid before calling this function.

Examples

using AlgebraicMultigrid

A = assemble(a; dirichlet = :boundary)
ml = amg_preconditioner(A)
P = aspreconditioner(ml)
sol = solve(LinearProblem(A, F), KrylovJL_CG(); Pl = P)

# Or, directly from the form:
uₕ = solve(a, l; dirichlet = bcs, preconditioner = :amg, solver = KrylovJL_CG())

See also linear_problem, assemble.

ilu_preconditioner(A::AbstractMatrix; kwargs...) -> ILUZero.ILU0Precon
ilu_preconditioner(a::Union{BilinearForm, Bramble.GhostBilinearForm}; dirichlet = nothing,
                    dirichlet_components = nothing, kwargs...) -> ILUZero.ILU0Precon

Build a zero-fill incomplete LU (ILU(0)) preconditioner for the matrix A – or for a assembled with the given Dirichlet conditions – ready to use directly as LinearSolve's Pl.

Convection-dominated forms (advection large relative to diffusion) assemble unsymmetric matrices far from an M-matrix, where algebraic multigrid (amg_preconditioner) does not degrade gracefully: measured on a 2D convection-diffusion system with diffusion 1e-2 against unit advection, ruge_stuben AMG failed to converge in 2000 GMRES iterations while ILU(0) converged in 18 (gpena/Bramble.jl#244). ILU(0) reuses A's own sparsity pattern for its factors, so it has no fill-in parameter to tune and is cheap to build, at the cost of a weaker preconditioner than a tuned incomplete factorization on harder systems.

Keywords

Forwarded to ILUZero.ilu0, which currently takes none beyond A itself.

The BilinearForm method matches assemble(a::BilinearForm; ...)'s own default and does not symmetrize; ILU(0) does not assume symmetry, so this is the natural entry point for the unsymmetric, convection-dominated forms this preconditioner targets.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

A form on ghost spaces (GhostGridSpace) is assembled to its reduced matrix, assemble(a; dirichlet, dirichlet_components), and factorised as any other.

Requires ILUZero.jl; call using ILUZero before calling this function.

Examples

using ILUZero

A = assemble(a; dirichlet = :boundary)
P = ilu_preconditioner(A)
sol = solve(LinearProblem(A, F), KrylovJL_GMRES(); Pl = P)

# Or, directly from the form:
uₕ = solve(a, l; dirichlet = bcs, preconditioner = :ilu0, solver = KrylovJL_GMRES())

See also amg_preconditioner, linear_problem, assemble.

AbstractMatrixFreePreconditioner{T}

A preconditioner P ≈ A for the matrix A of a BilinearForm, applied as P⁻¹ without storing A. A subtype supplies ldiv!(y, P, x) (y = P⁻¹ x) and ldiv!(P, x) (in place); P \ x returns P⁻¹ x in a new vector. That is the interface LinearSolve's Pl and Pr keywords take, so any subtype preconditions a Krylov solve of assemble(a) or of matrix_free_operator(a).

Type parameters

• T: The element type of P⁻¹ x for a vector x of element type T.

See also: JacobiPreconditioner, jacobi_preconditioner.

jacobi_preconditioner(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy = execution_policy(trial_space(a))) -> JacobiPreconditioner
jacobi_preconditioner(a::Bramble.GhostBilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy) -> JacobiPreconditioner
jacobi_preconditioner(op::MatrixFreeOperator) -> JacobiPreconditioner
jacobi_preconditioner(op::Bramble.GhostMatrixFreeOperator) -> JacobiPreconditioner

The Jacobi preconditioner of A = assemble(a; dirichlet, dirichlet_components), or of the matrix op stands for, from its diagonal alone: one walk of the form's stencil keeps the entries with row == col, and A is never built. Construction allocates the diagonal and the Dirichlet row mask, both of length ndofs. The walk runs under the operator's policy: under CpuThreaded or CpuPolyester its bands each own their rows, so the diagonal is the same, bit for bit, under every CPU policy at one thread count.

On a form over ghost spaces (GhostGridSpace), A is the reduced matrix E_vᵀ A_pad E_u that assemble(a) returns, and op the operator matrix_free_operator(a) returns for it: one walk of the padded form sums, for each degree of freedom, the padded entries whose row and column closures both reach it, so neither A_pad nor A is built. policy defaults to the padded space's. Under every policy, CpuSerial included, the walk keeps one vector of length ndofs per band of the threaded cut, which is why the diagonal is the same bits under each.

Arguments

• a: A square bilinear form (ndofs of trial and test space equal), on scalar, composite or ghost spaces.

• op: A square MatrixFreeOperator, or the operator matrix_free_operator returns for a ghost-space form; its Dirichlet rows are kept.

Keywords

• dirichlet, dirichlet_components, policy: As in matrix_free_operator.

Returns

• JacobiPreconditioner: ldiv!(y, P, x) gives x ./ diag(A).

Throws

• ArgumentError: policy is a GpuPolicy; device execution is tracked on milestone v4.7.0.

• ArgumentError: dirichlet_components names a leaf the test space does not have.

• DimensionMismatch: A is not square.

Examples

using Bramble, LinearAlgebra
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11, false))
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
P = jacobi_preconditioner(a; dirichlet = :boundary)
x = rand(ndofs(Wₕ))
P \ x ≈ x ./ diag(assemble(a; dirichlet = :boundary))

# output
true

See also: JacobiPreconditioner, matrix_free_operator.

JacobiPreconditioner{T, V <: AbstractVector{T}} <: AbstractMatrixFreePreconditioner{T}

The Jacobi preconditioner P = diag(A): ldiv!(y, P, x) sets y[i] = x[i] / A[i, i]. Stores the reciprocal diagonal, so an application is one multiply per entry and allocates nothing. A Dirichlet row of A is the identity row, so its entry is 1. Build one with jacobi_preconditioner.

A zero diagonal entry gives an infinite reciprocal, as x ./ diag(A) does.

Type parameters

• T: The element type of A.

• V: The vector type of the stored reciprocal diagonal.

See also: jacobi_preconditioner, AbstractMatrixFreePreconditioner.

chebyshev_preconditioner(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy = execution_policy(trial_space(a)), degree = 4, λmax = nothing, ratio = 30) -> ChebyshevPreconditioner
chebyshev_preconditioner(a::Bramble.GhostBilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy, degree = 4, λmax = nothing, ratio = 30) -> ChebyshevPreconditioner
chebyshev_preconditioner(op::MatrixFreeOperator; degree = 4, λmax = nothing, ratio = 30) -> ChebyshevPreconditioner
chebyshev_preconditioner(op::Bramble.GhostMatrixFreeOperator; degree = 4, λmax = nothing, ratio = 30) -> ChebyshevPreconditioner

The Chebyshev preconditioner of A = assemble(a; dirichlet, dirichlet_components), or of the matrix op stands for. It is a fixed polynomial in D⁻¹A, with D = diag(A), on the interval [λmax / ratio, λmax]. It is applied through products with the matrix-free operator, and D is read off one stencil walk as by jacobi_preconditioner, so A is never built. A must be symmetric positive definite, and every eigenvalue of D⁻¹A at most λmax: one above it can make the preconditioner indefinite. With dirichlet, use it in conjugate gradients only for a right-hand side that vanishes on the Dirichlet rows (see ChebyshevPreconditioner).

Eigenvalues below λmax / ratio are damped less, not lost, so ratio trades how far the polynomial reaches into the low spectrum against how flat it is on its interval. Construction allocates three vectors of length ndofs and, when λmax is not given, the two of max_eigenvalue_estimate.

On a form over ghost spaces (GhostGridSpace), A is the reduced matrix assemble(a) returns, applied by matrix_free_operator(a), and D its diagonal as jacobi_preconditioner reads it; policy defaults to the padded space's. That operator applies one product at a time, so neither it nor P may be shared between tasks.

Arguments

• a: A square bilinear form (ndofs of trial and test space equal), on scalar, composite or ghost spaces, whose matrix is SPD.

• op: A square MatrixFreeOperator, or the operator matrix_free_operator returns for a ghost-space form; its Dirichlet rows are kept.

Keywords

• dirichlet, dirichlet_components, policy: As in matrix_free_operator.

• degree: The number of Chebyshev steps, the degree of the residual polynomial; an application costs degree - 1 products with A (4 by default).

• λmax: The top of the interval, an upper bound for the spectrum of D⁻¹A, used as given (default: nothing, which takes max_eigenvalue_estimate(op; preconditioner = J) with J the Jacobi preconditioner of op, carrying its own safety factor).

• ratio: The ratio of the top of the interval to its bottom (30 by default).

Returns

• ChebyshevPreconditioner: ldiv!(y, P, x) gives p(D⁻¹A) D⁻¹ x.

Throws

• ArgumentError: degree < 1, ratio <= 1, or λmax (given or estimated) not finite and positive.

• ArgumentError: policy is a GpuPolicy; device execution is tracked on milestone v4.7.0.

• ArgumentError: dirichlet_components names a leaf the test space does not have.

• DimensionMismatch: A is not square.

Examples

With degree = 1 the polynomial is the constant \(1 / θ\), so P is Jacobi scaled by \(2 / (λ_{\max} + λ_{\min})\).

using Bramble, LinearAlgebra
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11, false))
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
P = chebyshev_preconditioner(a; degree = 1, λmax = 2.0, ratio = 3)
x = rand(ndofs(Wₕ))
P \ x ≈ (x ./ diag(assemble(a))) ./ (4 / 3)

# output
true

See also: ChebyshevPreconditioner, max_eigenvalue_estimate, jacobi_preconditioner.

ChebyshevPreconditioner{T, Op, V <: AbstractVector{T}} <: AbstractMatrixFreePreconditioner{T}

The Chebyshev preconditioner of an SPD operator A, scaled by its diagonal D: ldiv!(y, P, x) sets y = p(D⁻¹A) D⁻¹ x, the result of degree steps of Jacobi-preconditioned Chebyshev iteration for A y = x from y = 0 on the interval \([λ_{\min}, λ_{\max}]\). The residual polynomial \(1 - λ p(λ)\) is the scaled Chebyshev polynomial \(T_k((θ - λ) / δ) / T_k(θ / δ)\), with \(θ = (λ_{\max} + λ_{\min}) / 2\), \(δ = (λ_{\max} - λ_{\min}) / 2\) and \(k\) the degree, whose modulus is below 1 on \((0, λ_{\max}]\). So P⁻¹ is a fixed SPD operator whenever the spectrum of D⁻¹A lies in \((0, λ_{\max}]\), and plain conjugate gradients stays valid with it. Build one with chebyshev_preconditioner.

With Dirichlet rows A is not symmetric (only its rows are replaced), and neither is P⁻¹: the Dirichlet part of x feeds the interior part of y. Conjugate gradients with P is then valid only for a right-hand side and starting guess that vanish on the Dirichlet rows; for any other it can report convergence on a wrong solution.

An application costs degree - 1 products with A and allocates nothing: the two work vectors are stored in P, so one P must not be applied from two tasks at once.

Type parameters

• T: The element type of A.

• Op: The operator type, a MatrixFreeOperator or the operator matrix_free_operator returns for a ghost-space form.

• V: The vector type of the reciprocal diagonal and the work vectors.

See also: chebyshev_preconditioner, max_eigenvalue_estimate, AbstractMatrixFreePreconditioner.

fdm_preconditioner(a::BilinearForm; dirichlet = nothing) -> FDMPreconditioner

A preconditioner for assemble(a; dirichlet) that applies the fast-diagonalisation inverse of a's Laplacian-like part, factorised once, without assembling a's matrix. The Laplacian-like part keeps the terms fdm_solve can solve: with one mass per axis, every term that differs from the masses on at most one axis. Every term differing on two or more axes (a mixed derivative such as innerₕ(D₋ₓ(D₋ᵧ(u)), v), or a coefficient varying along two axes) is left out of the preconditioner, though the Krylov solver still sees it in assemble(a). For a form fdm_solve accepts, the preconditioner is the exact inverse.

Use it for a mixed-derivative form whose cross term the diffusion dominates, so that what it leaves out is small beside what it inverts. The larger the cross term, the more work it leaves to the Krylov solver. A form fdm_solve accepts needs no Krylov solver at all.

dirichlet is nothing (the unconstrained system) or :boundary (homogeneous Dirichlet on the whole mesh boundary): there the boundary rows of assemble(a; dirichlet = :boundary) are identity rows, so the preconditioner is the identity on them and the interior factorisation inside. Construction costs O(n_d^3) per axis; an ldiv! on host vectors costs O(N Σ_d n_d) and allocates nothing. The preconditioner is host-only.

Requires Kronecker.jl; call using Kronecker before calling this function.

Throws

• ArgumentError saying fdm_preconditioner does not support the form, with the reason, for a form with no usable Laplacian-like part: a is not separable, it is posed on a composite space, a mass is not symmetric positive definite, an axis has no term of its own once the two-axis terms are left out (the form has no Laplacian-like part), that part is singular, or a complex scale or coefficient makes the Kronecker operator's element type complex (a zero imaginary part included). These are fdm_solve's refusals, a two-axis term apart.

• ArgumentError: dirichlet is neither nothing nor :boundary.

Examples

using Bramble, Kronecker, LinearSolve
using Bramble: D₋ₓ, D₋ᵧ
Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (33, 25), (false, false))
Wₕ = gridspace(Ωₕ)
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)) +
                           0.25 * innerₕ(D₋ₓ(D₋ᵧ(u)), v))
P = fdm_preconditioner(a)
A = assemble(a)
F = rand(ndofs(Wₕ))
sol = solve(LinearProblem(A, F), KrylovJL_GMRES(); Pl = P)

See also: FDMPreconditioner, fdm_solve, jacobi_preconditioner.

FDMPreconditioner{T, F} <: AbstractMatrixFreePreconditioner{T}

The fast-diagonalisation preconditioner of a separable BilinearForm: ldiv!(y, P, x) applies the exact inverse of the form's Laplacian-like part (every term that differs from one mass per axis on at most one axis), factorised once, and allocates nothing on host vectors. Under dirichlet = :boundary it is the identity on the boundary rows. Build one with fdm_preconditioner, which needs using Kronecker; on a form over ghost spaces it inverts the Laplacian-like part of the reduced Kronecker operator.

Type parameters

• T: The element type of the form's Kronecker operator.

• F: The type of the stored factorisation, the one BrambleKroneckerExt builds.

See also: fdm_preconditioner, AbstractMatrixFreePreconditioner.

max_eigenvalue_estimate(op::MatrixFreeOperator; iterations = 20, preconditioner = nothing) -> T
max_eigenvalue_estimate(op::Bramble.GhostMatrixFreeOperator; iterations = 20, preconditioner = nothing) -> T

An upper estimate of the largest eigenvalue of the SPD matrix A that op stands for, or of P⁻¹A for an SPD preconditioner P, by power iteration on the matrix-free operator: iterations products w = A v (then w = P⁻¹ w) from a fixed pseudo-random start, each normalising v = w / ‖w‖, return 1.1 ‖w‖ for the last w.

For a symmetric A, ‖A v‖ ≤ λ_max for every unit v, so the bare power iterate is a lower bound, and it approaches λ_max slowly when the top eigenvalues cluster, as they do for a discrete diffusion operator. The factor 1.1 lifts it above λ_max once power iteration has come within 9%. On variable-coefficient diffusion on non-uniform meshes in one to three dimensions (up to 65² and 17³ points, five meshes each) the default iterations came within 8%, scaled by Jacobi or not, so the estimate landed between 1.02 λ_max and 1.1 λ_max. Overshooting is the side a Chebyshev preconditioner or smoother must err on: an eigenvalue above its interval can make it indefinite, while an interval a tenth too wide only shifts the polynomial slightly. With a preconditioner the operator P⁻¹A is not symmetric, only similar to a symmetric one, so ‖P⁻¹A v‖ may pass its top eigenvalue and the estimate overshoot a little more. Allocates two vectors of length ndofs.

Arguments

• op: A square MatrixFreeOperator of an SPD matrix, or the operator matrix_free_operator returns for a ghost-space form; its Dirichlet rows are kept.

Keywords

• iterations: The number of products with A (default: 20).

• preconditioner: nothing, or an AbstractMatrixFreePreconditioner P applied in place after each product, to estimate the top of the spectrum of P⁻¹A (default: nothing).

Returns

• T: The estimate, in the element type of op.

Throws

• ArgumentError: iterations < 1.

• DimensionMismatch: op is not square.

Examples

using Bramble, LinearAlgebra, Random
Random.seed!(1)
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 33, false))
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
λ = eigmax(Symmetric(Matrix(assemble(a))))
λ <= Bramble.max_eigenvalue_estimate(matrix_free_operator(a)) < 1.2λ

# output
true

See also: chebyshev_preconditioner, matrix_free_operator.

GeometricMeshHierarchy{D, M <: AbstractMeshType{D}}

The nested meshes of geometric multigrid, ordered coarse to fine: H[1] is the coarsest level and H[end] the mesh the hierarchy was built from. Each level holds every other point of the next finer one along every axis, and its markers are those of the finer level at the points it keeps. Every level has the type, and so the backend and element type, of the finest. Build one with GeometricMeshHierarchy(Ωₕ, levels).

A hierarchy indexes and iterates like a vector of meshes: length, getindex, keys, firstindex, lastindex and iteration are defined. It holds the finest mesh itself, not a copy, so it must be rebuilt after that mesh's points change (change_points!, iterative_refinement!): the coarser levels would no longer nest in it.

Type parameters

• D: The dimension of the meshes.

• M: The mesh type of every level.

See also: mesh, change_points!.

prolongate!(xf, H::GeometricMeshHierarchy, l::Integer, xc) -> xf

The prolongation xf = P xc from level l - 1 of H to level l: the multilinear interpolant of the coarse grid function xc at the points of H[l]. A fine point that is a coarse point takes its value; one between two coarse points \(x_{j}, x_{j+1}\) along an axis takes \(((x_{j+1} - x) u_j + (x - x_{j}) u_{j+1}) / (x_{j+1} - x_{j})\), per axis, from the actual (non-uniform) coordinates. P equals Bramble.interpolation_matrix(gridspace(H[l]), gridspace(H[l - 1])), which is never built.

The product runs as one sweep per axis, in place in xf, under the execution policy of the backend of H[l]. It allocates nothing on a CpuSerial backend; on CpuThreaded only the task spawns, whatever the grid size. It keeps no state, so transfers on one hierarchy may run from several tasks at once.

Arguments

• xf: The fine grid function, overwritten: a vector of length npoints(H[l]) or a VectorElement of a space on H[l].

• H: The hierarchy.

• l: The index of the fine level, 2 ≤ l ≤ length(H).

• xc: The coarse grid function: a vector of length npoints(H[l - 1]) or a VectorElement of a space on H[l - 1]. It must not alias xf.

Returns

• xf, holding P xc.

Throws

• ArgumentError: l is not in 2:length(H).

• DimensionMismatch: xf or xc has the wrong length.

• ArgumentError: xf or xc is not 1-based, or the two may alias.

• ArgumentError: the backend's policy is a GpuPolicy; device transfers are tracked on milestone v4.7.0.

Examples

Interpolation reproduces a linear function on a non-uniform mesh.

using Bramble
H = GeometricMeshHierarchy(mesh(domain(interval(0.0, 1.0)), 9, false), 2)
xf = zeros(npoints(H[2]))
prolongate!(xf, H, 2, 2 .* points(H[1]) .+ 1)
xf ≈ 2 .* points(H[2]) .+ 1

# output
true

See also: coarsen!, GeometricMeshHierarchy.

coarsen!(xc, H::GeometricMeshHierarchy, l::Integer, xf) -> xc

The transpose of prolongate!: xc = Pᵀ xf, from level l of H to level l - 1, with P the prolongation from level l - 1 to level l. There is no scaling: a coarse value is the sum of the fine values around it, each weighted by what the coarse point contributes to it under interpolation. Bramble's bilinear forms carry the discrete measure, so Pᵀ A_f P has the scaling of the coarse form's matrix, not that matrix itself, and they are equal for the 1D stiffness, even on a non-uniform mesh, but in 2D and 3D, or for a mass form, they differ entrywise by a quarter to a half. Full weighting (Pᵀ / 2^D), the choice for strong-form difference equations, would be off by 2^D.

Each coarse value is gathered from its \(3^D\) nearest fine points, with the weights prolongate! uses, under the execution policy of the backend of H[l]. It allocates nothing on a CpuSerial backend; on CpuThreaded only the task spawns, whatever the grid size. It keeps no state, so transfers on one hierarchy may run from several tasks at once.

Arguments

• xc: The coarse grid function, overwritten: a vector of length npoints(H[l - 1]) or a VectorElement of a space on H[l - 1].

• H: The hierarchy.

• l: The index of the fine level, 2 ≤ l ≤ length(H).

• xf: The fine grid function, not modified: a vector of length npoints(H[l]) or a VectorElement of a space on H[l]. It must not alias xc.

Returns

• xc, holding Pᵀ xf.

Throws

• ArgumentError: l is not in 2:length(H).

• DimensionMismatch: xf or xc has the wrong length.

• ArgumentError: xf or xc is not 1-based, or the two may alias.

• ArgumentError: the backend's policy is a GpuPolicy; device transfers are tracked on milestone v4.7.0.

Examples

The rows of P sum to one, so Pᵀ preserves the sum of a fine grid function.

using Bramble
H = GeometricMeshHierarchy(mesh(domain(interval(0.0, 1.0) × interval(0.0, 2.0)), (9, 5), false), 2)
xf = rand(npoints(H[2]))
xc = coarsen!(zeros(npoints(H[1])), H, 2, xf)
(length(xc), sum(xc) ≈ sum(xf))

# output
(15, true)

See also: prolongate!, GeometricMeshHierarchy.

AbstractSmoother{T}

A multigrid smoother for the matrix A of a BilinearForm: smooth!(s, x, b) moves the iterate x towards the solution of A x = b in place, damping the high-frequency part of its error. A subtype stores its operator, the Jacobi diagonal and its work vectors, so smooth! allocates nothing on a serial policy, and one smoother must not be used from two tasks at once.

All three smoothers are point smoothers: each unknown is updated from its own equation. They damp the high frequencies of an isotropic problem, but where cells are stretched the error that oscillates only along the long side of a cell is smooth for them and oscillatory for the coarse grid, and neither removes it. With one application before and after an exact Galerkin coarse correction, on mass plus variable diffusion, the two-grid contraction factors quoted for each smoother (uniform meshes) rose to 0.74–0.93 on a 2D mesh of 33² graded 7.4-fold along each axis (cell aspect ratio 6.9), and to 0.68–0.995 on the random non-uniform meshes of mesh(…, false) (2D 33², 3D 9³, six seeds), tracking the largest cell aspect ratio: 0.69 at 36, 0.995 at 27400. This held for every smoother and degree; in 1D, where no cell is stretched, they did not change.

Type parameters

• T: The element type of A.

See also: JacobiSmoother, ChebyshevSmoother, RedBlackGaussSeidel, smooth!.

JacobiSmoother{T, Op, V <: AbstractVector{T}} <: AbstractSmoother{T}

Damped Jacobi: one sweep of smooth! sets x ← x + ω D⁻¹ (b - A x), with D the diagonal of A. A sweep costs one product with A and allocates nothing: the residual vector is stored in the smoother, so one smoother must not be used from two tasks at once. Build one with jacobi_smoother.

Type parameters

• T: The element type of A.

• Op: The operator type, a MatrixFreeOperator or a ghost-space operator.

• V: The vector type of the reciprocal diagonal and the residual.

See also: jacobi_smoother, AbstractSmoother.

ChebyshevSmoother{T, Op, V <: AbstractVector{T}} <: AbstractSmoother{T}

Chebyshev smoothing: smooth! runs degree steps of Jacobi-scaled Chebyshev iteration from the iterate passed in, on the interval \([λ_{\max} / 4, λ_{\max}]\) of the spectrum of D⁻¹A (the recurrence of ChebyshevPreconditioner). The error x - A⁻¹ b is multiplied by \(T_k((θ - D⁻¹A) / δ) / T_k(θ / δ)\), with \(θ = 5 λ_{\max} / 8\), \(δ = 3 λ_{\max} / 8\) and \(k\) the degree, whose modulus on the interval is at most \(1 / T_k(5 / 3)\). An application costs degree products with A and allocates nothing: the two work vectors are stored in the smoother, so one smoother must not be used from two tasks at once. Build one with chebyshev_smoother.

Type parameters

• T: The element type of A.

• Op: The operator type, a MatrixFreeOperator or a ghost-space operator.

• V: The vector type of the reciprocal diagonal and the work vectors.

See also: chebyshev_smoother, AbstractSmoother.

RedBlackGaussSeidel{T, Op, V <: AbstractVector{T}, D} <: AbstractSmoother{T}

Red-black Gauss-Seidel on a scalar grid space: the grid points are coloured by the parity of the sum of their (zero-based) indices, red when it is even, so the first point is red. One sweep of smooth! updates the red points, then the black ones, each by x[i] ← x[i] + (b - A x)[i] / A[i, i] with the residual of the current iterate. As no entry of A couples two distinct points of one colour (checked by red_black_gauss_seidel), each half-sweep is exact Gauss-Seidel on its colour, and a sweep equals forward Gauss-Seidel in the red-then-black ordering. A sweep costs two products with A and allocates nothing: the residual vector is stored in the smoother, so one smoother must not be used from two tasks at once. smooth!(s, x, b; reverse = true) updates the black points first, and that sweep is backward Gauss-Seidel in the red-then-black ordering, the adjoint of the forward one, and is what a symmetric multigrid cycle post-smooths with. The two-grid factors quoted below are for red then black on both sides.

Type parameters

• T: The element type of A.

• Op: The operator type, a MatrixFreeOperator.

• V: The vector type of the reciprocal diagonal and the residual.

• D: The dimension of the grid.

See also: red_black_gauss_seidel, AbstractSmoother.

jacobi_smoother(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy = execution_policy(trial_space(a)), ω = nothing, sweeps = 1) -> JacobiSmoother
jacobi_smoother(a::Bramble.GhostBilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy, ω = nothing, sweeps = 1) -> JacobiSmoother
jacobi_smoother(op::MatrixFreeOperator; ω = nothing, sweeps = 1) -> JacobiSmoother
jacobi_smoother(op::Bramble.GhostMatrixFreeOperator; ω = nothing, sweeps = 1) -> JacobiSmoother

The damped Jacobi smoother of A = assemble(a; dirichlet, dirichlet_components), or of the matrix op stands for: each sweep sets x ← x + ω D⁻¹ (b - A x), with D = diag(A) read off one stencil walk as by jacobi_preconditioner. Construction allocates two vectors of length ndofs and no product with A, so it is cheap enough to build once per multigrid level. A Dirichlet row of A is the identity row, so a sweep moves x[i] a fraction ω of the way to b[i] there.

For a diagonally dominant form, such as second-order diffusion plus mass, the spectrum of D⁻¹A lies in (0, 2], and a sweep multiplies an eigencomponent of the error by 1 - ωλ. In d dimensions the components a grid coarsened by 2 cannot represent have λ in about [1/d, 2] (the 3-, 5- and 7-point Laplacians), and the default ω = 2d / (2d + 1) (2/3, 4/5, 6/7) is the weight that minimises |1 - ωλ| over that band. With one sweep before and after an exact Galerkin coarse correction, on mass plus variable diffusion, it gave two-grid contraction factors 0.36 (2D, 33²) and 0.51 (3D, 9³), against 0.44 and 0.56 for ω = 2/3; in 1D the two are the same weight (0.11 on 65 points, uniform, graded or random). Mesh non-uniformity with large cell aspect ratios slows every point smoother (see AbstractSmoother).

Arguments

• a: A square bilinear form (ndofs of trial and test space equal), on scalar, composite or ghost spaces; on ghost spaces policy defaults to the padded space's.

• op: A square MatrixFreeOperator, or the operator matrix_free_operator returns for a form on ghost spaces, whose diagonal is that of the reduced matrix assemble(a); its Dirichlet rows are kept.

Keywords

• dirichlet, dirichlet_components, policy: As in matrix_free_operator.

• ω: The damping weight (default: nothing, for 2d / (2d + 1) with d the dimension of the test space).

• sweeps: The number of sweeps per smooth! (default: 1, so that a multigrid cycle's pre- and post-smoothing counts are sweep counts).

Returns

• JacobiSmoother.

Throws

• ArgumentError: sweeps < 1, or ω not finite and positive.

• ArgumentError: policy is a GpuPolicy; device execution is tracked on milestone v4.7.0.

• ArgumentError: dirichlet_components names a leaf the test space does not have.

• DimensionMismatch: A is not square.

Examples

In 1D the default weight is 2/3.

using Bramble, LinearAlgebra
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11, false))
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
A = assemble(a)
x, b = rand(ndofs(Wₕ)), rand(ndofs(Wₕ))
x₀ = copy(x)
smooth!(jacobi_smoother(a), x, b) ≈ x₀ .+ (2 / 3) .* (b .- A * x₀) ./ diag(A)

# output
true

See also: JacobiSmoother, smooth!, chebyshev_smoother, red_black_gauss_seidel.

chebyshev_smoother(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy = execution_policy(trial_space(a)), degree = 2, λmax = nothing) -> ChebyshevSmoother
chebyshev_smoother(a::Bramble.GhostBilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy, degree = 2, λmax = nothing) -> ChebyshevSmoother
chebyshev_smoother(op::MatrixFreeOperator; degree = 2, λmax = nothing) -> ChebyshevSmoother
chebyshev_smoother(op::Bramble.GhostMatrixFreeOperator; degree = 2, λmax = nothing) -> ChebyshevSmoother

The Chebyshev smoother of A = assemble(a; dirichlet, dirichlet_components), or of the matrix op stands for: smooth! runs degree steps of the Jacobi-scaled Chebyshev recurrence of chebyshev_preconditioner from the iterate passed in, on the interval [λmax / 4, λmax] of the spectrum of D⁻¹A, D = diag(A). The interval is the upper part of the spectrum a coarser level cannot represent, so the smoother damps all of it and leaves the rest to the coarse correction. A must be SPD with every eigenvalue of D⁻¹A at most λmax; one above the interval is amplified.

Each step damps the upper spectrum by more than a Jacobi sweep does for the same product: on the interval the error contracts by at most 1 / T_k(5/3) in k products, 0.22 at degree 2. With one application before and after an exact Galerkin coarse correction, on mass plus variable diffusion on uniform meshes, degree 2 gave two-grid contraction factors 0.055 (1D, 65 points), 0.088 (2D, 33²) and 0.20 (3D, 9³) for four products with A per cycle, against 0.062, 0.20 and 0.33 for two sweeps of Jacobi with ω = 2/3 at the same cost; degree 3 gave 0.034, 0.064 and 0.098 for six. Per product, degree 2 is the fastest or tied in every dimension, and on a graded 2D mesh too (0.82, against 0.74 at degree 3), so it is the default.

Construction allocates four vectors of length ndofs and, when λmax is not given, runs max_eigenvalue_estimate's 20 products with A and allocates its two vectors.

Arguments

• a: A square bilinear form (ndofs of trial and test space equal), on scalar, composite or ghost spaces, whose matrix is SPD; on ghost spaces policy defaults to the padded space's.

• op: A square MatrixFreeOperator, or the operator matrix_free_operator returns for a form on ghost spaces, whose diagonal is that of the reduced matrix assemble(a); its Dirichlet rows are kept.

Keywords

• dirichlet, dirichlet_components, policy: As in matrix_free_operator.

• degree: The number of Chebyshev steps, and of products with A, per smooth! (default: 2).

• λmax: The top of the interval, an upper bound for the spectrum of D⁻¹A, used as given (default: nothing, for max_eigenvalue_estimate(op; preconditioner = J) with J the Jacobi preconditioner of op).

Returns

• ChebyshevSmoother.

Throws

• ArgumentError: degree < 1, or λmax (given or estimated) not finite and positive.

• ArgumentError: policy is a GpuPolicy; device execution is tracked on milestone v4.7.0.

• ArgumentError: dirichlet_components names a leaf the test space does not have.

• DimensionMismatch: A is not square.

Examples

At degree 1 the smoother is damped Jacobi with ω = 1 / θ = 8 / (5 λmax).

using Bramble
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11, false))
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
x, b = rand(ndofs(Wₕ)), rand(ndofs(Wₕ))
y = copy(x)
smooth!(chebyshev_smoother(a; degree = 1, λmax = 2.0), x, b) ≈ smooth!(jacobi_smoother(a; ω = 0.8), y, b)

# output
true

See also: ChebyshevSmoother, smooth!, chebyshev_preconditioner.

red_black_gauss_seidel(a::BilinearForm; dirichlet = nothing, policy = execution_policy(trial_space(a)), sweeps = 1) -> RedBlackGaussSeidel
red_black_gauss_seidel(a::Bramble.GhostBilinearForm; dirichlet = nothing, policy, sweeps = 1) -> RedBlackGaussSeidel
red_black_gauss_seidel(op::MatrixFreeOperator; sweeps = 1) -> RedBlackGaussSeidel

The red-black Gauss-Seidel smoother of A = assemble(a; dirichlet), or of the matrix op stands for (see RedBlackGaussSeidel): each sweep updates the points whose zero-based indices sum to an even number, then the rest, by a masked Jacobi step with the current residual, which is exact Gauss-Seidel on each colour because no entry of A couples two distinct points of one colour. That holds for the 3-, 5- and 7-point forms of second-order diffusion and for mass terms; a form with a mixed difference such as D₋ᵧ(D₋ₓ(u)) (a 9-point box) couples diagonal neighbours, which share a colour, and is refused. The check reads the offsets of every term's trial and test sides (stencil_offsets) and ignores collapsed axes (one point), where no entry lands off the axis. A Dirichlet row of A is the identity row, so a sweep sets x[i] = b[i] there.

Construction allocates two vectors of length ndofs and the small offset lists, and no product with A, so it is cheap enough to build once per multigrid level. Red-black is defined on one grid: a composite space, whose leaves couple at a shared point, is refused. With one sweep before and after an exact Galerkin coarse correction, on mass plus variable diffusion on uniform meshes, it gave two-grid contraction factors 0.002 (1D, 65 points), 0.063 (2D, 33²) and 0.18 (3D, 9³) for four products with A per cycle, close to degree-2 chebyshev_smoother at the same cost. The default sweeps = 1 makes a multigrid cycle's pre- and post-smoothing counts sweep counts.

Arguments

• a: A square bilinear form on a scalar grid space, with the same mesh for trial and test space and no interpolation (πₕ). A form on ghost spaces is refused, as its operator is.

• op: A MatrixFreeOperator of such a form; its Dirichlet rows are kept. The operator of a form on ghost spaces is refused.

Keywords

• dirichlet, policy: As in matrix_free_operator.

• sweeps: The number of red-then-black sweeps per smooth! (default: 1).

Returns

• RedBlackGaussSeidel.

Throws

• ArgumentError: sweeps < 1.

• ArgumentError: the form couples two distinct points of one colour, is on a composite space, has an interpolation, or its trial and test meshes differ in shape.

• ArgumentError: a is a form on ghost spaces, or op its operator, whose closure couples a facet node to nodes of its own colour; use jacobi_smoother or chebyshev_smoother.

• ArgumentError: policy is a GpuPolicy; device execution is tracked on milestone v4.7.0.

• DimensionMismatch: A is not square.

Examples

After a sweep the residual vanishes at the black points, whose equations the black half has just solved exactly.

using Bramble, Random
Random.seed!(1)
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11, false))
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
x, b = rand(ndofs(Wₕ)), rand(ndofs(Wₕ))
smooth!(red_black_gauss_seidel(a), x, b)
maximum(abs, (b - assemble(a) * x)[2:2:end]) < 1e-12

# output
true

See also: RedBlackGaussSeidel, smooth!, jacobi_smoother.

smooth!(s::AbstractSmoother, x, b) -> x
smooth!(s::RedBlackGaussSeidel, x, b; reverse::Bool = false) -> x

Apply the smoother s to the iterate x for A x = b, in place: the sweeps of a JacobiSmoother or RedBlackGaussSeidel, the Chebyshev steps of a ChebyshevSmoother. Allocates nothing on a CpuSerial policy; on a threaded one, only the task launches of each product with A. The work vectors live in s, so one smoother must not be used from two tasks at once.

The Jacobi and Chebyshev error propagators are polynomials in D⁻¹A, self-adjoint in the A inner product. A red-black sweep is forward Gauss-Seidel in the red-then-black ordering, and reverse = true runs each sweep black then red, its adjoint (backward Gauss-Seidel in the same ordering).

Arguments

• s: The smoother.

• x: The iterate, overwritten: a vector of length ndofs, or a VectorElement.

• b: The right-hand side, not modified: a vector of length ndofs, or a VectorElement. It must not alias x.

Keywords

• reverse: Red-black only: update the black points before the red ones in every sweep (default: false).

Returns

• x, the smoothed iterate.

Throws

• DimensionMismatch: x or b has the wrong length.

• ArgumentError: x or b is not 1-based, or the two may alias.

Examples

One damped Jacobi sweep with ω = 1/2 on the 1D Laplacian with its highest-frequency error (-1)^i: away from the ends D⁻¹A doubles it, so the sweep removes it there exactly.

using Bramble
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11, true))
a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))
x = [(-1.0)^i for i in 1:11]
smooth!(jacobi_smoother(a; ω = 0.5), x, zeros(11))
maximum(abs, x[2:10]) < 1e-12

# output
true

See also: AbstractSmoother, jacobi_smoother, chebyshev_smoother, red_black_gauss_seidel.

gmg_preconditioner(build, Ωₕ::AbstractMeshType; levels = nothing, cycle = :V, ν₁ = 2, ν₂ = 2, smoother = op -> chebyshev_smoother(op)) -> GMGPreconditioner

The geometric multigrid preconditioner of the form build(gridspace(Ωₕ)). ldiv!(y, P, x) runs one cycle for A y = x from y = 0 over the hierarchy GeometricMeshHierarchy(Ωₕ, levels). Each level l is rediscretised, so its form is build(gridspace(H[l])), applied through matrix_free_operator under the backend's execution policy, and smoothed by smoother(op). Levels are joined by prolongate! and coarsen! (Pᵀ, unscaled, since Bramble's forms carry the discrete measure). The coarsest level's matrix is assembled, stored dense and LU-factorised once, and solved exactly.

A V-cycle on level l smooths ν₁ times, restricts the residual, runs one cycle on level l - 1 from zero (two for a W-cycle), adds the prolongated correction and smooths ν₂ times with the adjoint smoother (a RedBlackGaussSeidel runs black then red). With ν₁ == ν₂ and an SPD A, the V- and W-cycles are symmetric positive definite, and conjugate gradients may use them (Pl = P in LinearSolve's KrylovJL_CG). The full multigrid cycle :FMG (fmg!) is a linear operator but not symmetric.

build(W) must return a square BilinearForm on the scalar grid space W it is given, with no Dirichlet rows (a form carries none; A is assemble(build(W))), or a square form on ghost spaces built on mesh(W), such as form(Cₕ, Cₕ, f) with Cₕ = constrained(gridspace(mesh(W); ghosts = ...), relations). On ghost spaces, every level has ghosts on the same facets and every ghost is eliminated by a relation; Dirichlet relations u => g are allowed, and the unknowns of each level are its ndofs. Neumann, Robin, mixed and pure Neumann relations all work: when the constant lies in the kernel of the coarsest matrix, its solve is bordered by the constant and returns the mean-zero solution, so a compatible right-hand side (b = A x) converges, to some solution of A u = b. The smoothers read the reduced diagonal, and red-black Gauss-Seidel is refused. Conjugate gradients needs A SPD, such as mass plus diffusion with natural boundary conditions. A Dirichlet problem needs its constrained rows eliminated or lifted into the right-hand side before a form without them is passed here. Every grid function in the form (a coefficient Rₕ(W, κ), say) must be built from W inside build: one built on the finest space and captured is accepted on every level and gives a wrong coarse operator.

On meshes whose cells have bounded aspect ratio, CG iteration counts do not grow with the mesh: with the defaults, mass plus variable diffusion took 6 iterations (relative residual 1e-8) at every size from 2D 33² to 513², and 7 from 3D 17³ to 129³, on uniform points each jittered by up to ±0.3h; the W-cycle took 5 and full multigrid 4 at 65². The smoothers are point smoothers, and they stall where cells are stretched (see AbstractSmoother). On the random meshes of mesh(…, false), whose largest cell aspect ratio grows with n (96 at 33², 52600 at 513²), CG took 14–25, 26–37 and 32–116 iterations at 2D 33², 65² and 129² over four draws, and 87 at 513² in one. A random base mesh refined with iterative_refinement! keeps its aspect ratio, but each stretched base cell becomes a block of stretched cells, and iterations still grow per level: 11, 16, 21, 26, 26 from 17² to 257² from a 2D base of 9² (aspect ratio 10.9); 23 to 56 over 33² to 257² from a base of 17² (aspect ratio 43); 14 to 29 over 9³ to 65³ from a 3D base of 5³ (aspect ratio 15); 21 to 41 over 17³ to 65³ from a base of 9³ (aspect ratio 36). Red-black Gauss-Seidel gave the same counts.

Construction builds one form, operator and smoother per level (the default smoother runs max_eigenvalue_estimate's 20 products per level), assembles the coarsest matrix, and allocates three work vectors per level. The work vectors live in P, so one preconditioner must not be applied from two tasks at once. On a threaded policy an application allocates only the task launches of its products with A and its transfers: that grows with the number of products (a W-cycle doubles the cycles per level) but not with the grid size.

Arguments

• build: A function of a ScalarGridSpace W returning a square BilinearForm on W, or a square form on ghost spaces built on mesh(W).

• Ωₕ: The finest mesh.

Keywords

• levels: The number of levels, Ωₕ included, an integer (nothing by default, for as many as coarsening by 2 allows while every axis with more than one point keeps at least three, so 33² gives 5 levels down to 3², and 97 points, with 96 = 3 ⋅ 2⁵, give 6 down to 4).

• cycle: :V, :W or :FMG, the cycle ldiv! runs (:V by default).

• ν₁, ν₂: The number of smooth! calls before and after the coarse correction (2 each by default).

• smoother: A function of a level's MatrixFreeOperator returning its AbstractSmoother (default: op -> chebyshev_smoother(op), degree 2, which damped most per product with A in the two-grid measurements of chebyshev_smoother). A smoother other than the three Bramble provides is post-smoothed with smooth! as well, so it must be self-adjoint in the A inner product for the cycle to be symmetric.

Returns

• GMGPreconditioner.

Throws

• ArgumentError: cycle is not :V, :W or :FMG, ν₁ or ν₂ is negative, or ν₁ + ν₂ == 0 (no smoothing leaves the cycle singular).

• ArgumentError: levels is not an integer, or is invalid for Ωₕ (see GeometricMeshHierarchy).

• ArgumentError: build(W) does not return a square form on W or on ghost spaces on mesh(W), returns a form on ghost spaces on some levels only, or with ghosts on other facets than on the finest level.

• ArgumentError: on ghost spaces, a ghost is kept as an unknown (a ghosted facet without a relation), or the trial and test spaces have different degrees of freedom.

• ArgumentError: build(W) returns a form over polygon spaces (polygon_spaces), and so does gmg_preconditioner(a) on one: multigrid on them is gpena/Bramble.jl#846.

• ArgumentError: smoother(op) does not return an AbstractSmoother.

• ArgumentError: the coarsest level has more than 4096 points: Ωₕ's point counts minus one need a larger power of two as a factor.

• ArgumentError: the backend's policy is a GpuPolicy; device execution is tracked on milestone v4.7.0.

Examples

The V-cycle is symmetric.

using Bramble, LinearAlgebra
Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (33, 33), (true, true))
P = gmg_preconditioner(W -> form(W, W, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v))), Ωₕ)
x, y = rand(33^2), rand(33^2)
(length(P.hierarchy), dot(P \ x, y) ≈ dot(x, P \ y))

# output
(5, true)

See also: GMGPreconditioner, gmg_solve, v_cycle!, GeometricMeshHierarchy.

GMGPreconditioner{T, H, Op, S, F, V <: AbstractVector{T}, G} <: AbstractMatrixFreePreconditioner{T}

A geometric multigrid cycle as a preconditioner: ldiv!(y, P, x) sets y to one cycle for A y = x from y = 0, where A is the finest level's form. It holds the GeometricMeshHierarchy, one MatrixFreeOperator and one smoother per level, the dense factorisation of the coarsest level's matrix, and the per-level work vectors (of each level's ndofs), so an application allocates nothing on a CpuSerial policy and one P must not be applied from two tasks at once. Build one with gmg_preconditioner; the cycles themselves are v_cycle!, w_cycle! and fmg!.

Type parameters

• T: The element type of A.

• H: The hierarchy type.

• Op: The operator type of the levels.

• S: The smoother type of the levels.

• F: The type of the coarsest level's factorisation.

• V: The vector type of the work vectors.

• G: Nothing on plain spaces; on ghost spaces, the type of the map from each level's degrees of freedom to its physical nodes, which the transfers go through.

See also: gmg_preconditioner, gmg_solve, AbstractMatrixFreePreconditioner.

gmg_solve(build, Ωₕ::AbstractMeshType, b; levels = nothing, cycle = :V, ν₁ = 2, ν₂ = 2, smoother = op -> chebyshev_smoother(op), tol = nothing, maxiters = 100) -> VectorElement

The solution of A u = b by multigrid iteration, with A the matrix of build(gridspace(Ωₕ)): from u = 0, cycles of gmg_preconditioner(build, Ωₕ; levels, cycle, ν₁, ν₂, smoother) repeat until the residual satisfies ‖b - A u‖ ≤ tol ‖b‖ (2-norm). A :V or :W cycle is repeated as itself; :FMG runs fmg! once and V-cycles after it. Each cycle costs one more product with A for the residual.

This is a stationary iteration: its contraction per cycle is the cycle's, so where point smoothers stall (stretched cells, see gmg_preconditioner), conjugate gradients preconditioned by the V-cycle converges in fewer cycles. A needs no symmetry here, but the cycle must contract. As in gmg_preconditioner, every grid function in the form must be built from W inside build.

The residual floors at a multiple of the rounding error of the element type T of A: in Float32 it stalled near 1e-5 relative at 65², so tol must stay above that. The default √eps(T) (1.5e-8 in Float64, 3.5e-4 in Float32) does.

Arguments

• build: A function of a ScalarGridSpace W returning a square BilinearForm on W, or one on ghost spaces built on mesh(W), as in gmg_preconditioner.

• Ωₕ: The finest mesh.

• b: The right-hand side: a vector of length ndofs of the finest form's spaces (npoints(Ωₕ) on a plain space) or a VectorElement, such as assemble(l) for a linear form l.

Keywords

• levels, cycle, ν₁, ν₂, smoother: As in gmg_preconditioner.

• tol: The relative residual reached (default: nothing, for √eps(T) with T the element type of A).

• maxiters: The largest number of cycles (default: 100).

Returns

• VectorElement: The solution, an element of the trial space of the finest form.

Throws

• As gmg_preconditioner.

• DimensionMismatch: b has the wrong length.

• ArgumentError: b is not 1-based, tol is not positive, or maxiters < 1.

• ErrorException: the residual is above tol ‖b‖ after maxiters cycles; the message names tol and T, since a tol near the rounding floor of T is never reached.

Examples

using Bramble, LinearAlgebra
Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (33, 33), (true, true))
build(W) = form(W, W, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
b = rand(33^2)
uₕ = gmg_solve(build, Ωₕ, b; tol = 1e-10)
norm(assemble(build(gridspace(Ωₕ))) * parent(uₕ) - b) <= 1e-10 * norm(b)

# output
true

See also: gmg_preconditioner, v_cycle!, fmg!.

v_cycle!(x, P::GMGPreconditioner, b) -> x

One V-cycle for A x = b from the iterate x, in place, with A the finest level's form of P: ν₁ smoothing steps, the coarse correction by one V-cycle on the next level from zero (an exact solve on the coarsest), then ν₂ adjoint smoothing steps. Allocates nothing on a CpuSerial policy. The work vectors live in P, so it must not run on one P from two tasks at once.

Arguments

• x: The iterate, overwritten: a vector of length ndofs, or a VectorElement.

• P: The multigrid preconditioner.

• b: The right-hand side, not modified: a vector of length ndofs, or a VectorElement. It must not alias x.

Returns

• x, the iterate after the cycle.

Throws

• DimensionMismatch: x or b has the wrong length.

• ArgumentError: x or b is not 1-based, or the two may alias.

Examples

Each V-cycle divides the residual by more than five.

using Bramble, LinearAlgebra
using Bramble: v_cycle!
Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (33, 33), (true, true))
build(W) = form(W, W, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
P = gmg_preconditioner(build, Ωₕ)
A = assemble(build(gridspace(Ωₕ)))
b, x = rand(33^2), zeros(33^2)
r₀ = norm(b)
v_cycle!(x, P, b)
r₁ = norm(b - A * x)
v_cycle!(x, P, b)
(r₁ < r₀ / 5, norm(b - A * x) < r₁ / 5)

# output
(true, true)

See also: w_cycle!, fmg!, gmg_preconditioner.

w_cycle!(x, P::GMGPreconditioner, b) -> x

One W-cycle for A x = b from the iterate x, in place: as v_cycle!, but the coarse correction on every level below the finest runs two cycles in succession (the second from the first's result) instead of one. Allocates nothing on a CpuSerial policy.

Arguments

• x: The iterate, overwritten: a vector of length ndofs, or a VectorElement.

• P: The multigrid preconditioner.

• b: The right-hand side, not modified: a vector of length ndofs, or a VectorElement. It must not alias x.

Returns

• x, the iterate after the cycle.

Throws

• DimensionMismatch: x or b has the wrong length.

• ArgumentError: x or b is not 1-based, or the two may alias.

Examples

On two levels the coarse solve is exact, so the W- and V-cycles agree.

using Bramble
using Bramble: v_cycle!, w_cycle!
Ωₕ = mesh(domain(interval(0.0, 1.0)), 17, true)
P = gmg_preconditioner(W -> form(W, W, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v))), Ωₕ; levels = 2)
b = rand(17)
w_cycle!(zeros(17), P, b) ≈ v_cycle!(zeros(17), P, b)

# output
true

See also: v_cycle!, fmg!, gmg_preconditioner.

fmg!(x, P::GMGPreconditioner, b) -> x

Full multigrid for A x = b, overwriting x (its input is not used): b is restricted by coarsen! to every level, the coarsest is solved exactly, and on each finer level the prolongated solution of the level below starts one V-cycle. The result is linear in b but not symmetric. Allocates nothing on a CpuSerial policy.

Arguments

• x: The solution, overwritten: a vector of length ndofs, or a VectorElement.

• P: The multigrid preconditioner.

• b: The right-hand side, not modified: a vector of length ndofs, or a VectorElement. It must not alias x.

Returns

• x, the full multigrid approximation of A⁻¹ b.

Throws

• DimensionMismatch: x or b has the wrong length.

• ArgumentError: x or b is not 1-based, or the two may alias.

Examples

On a smooth right-hand side, one full multigrid cycle leaves a residual more than ten times smaller than one V-cycle from zero does.

using Bramble, LinearAlgebra
using Bramble: v_cycle!, fmg!
Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (33, 33), (true, true))
build(W) = form(W, W, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
P = gmg_preconditioner(build, Ωₕ)
A = assemble(build(gridspace(Ωₕ)))
b = A * vec([sin(3x) * cos(2y) for x in range(0, 1, 33), y in range(0, 1, 33)])
r_fmg = norm(b - A * fmg!(zeros(33^2), P, b))
r_v = norm(b - A * v_cycle!(zeros(33^2), P, b))
r_fmg < r_v / 10

# output
true

See also: v_cycle!, w_cycle!, gmg_preconditioner.

sparse_factorize(A::SparseMatrixCSC; solver::Symbol = :default, sym = :auto, kwargs...) -> Factorization
sparse_factorize(a::BilinearForm; solver::Symbol = :default, sym = :auto,
                 dirichlet = nothing, dirichlet_components = nothing, kwargs...) -> Factorization

Compute the sparse direct factorization of A (or the assembled matrix of a) using the requested solver backend.

Solvers (solver)

• :default or :suitesparse: SuiteSparse (CHOLMOD for SPD/symmetric, UMFPACK for unsymmetric).

• :accelerate: Apple Accelerate native libSparse on macOS (Cholesky, LDLᵀ, LUTPP, QR).

• :mumps: MUMPS multifrontal parallel direct solver.

• :sparspak: pure-Julia sparse direct LU, zero binary dependencies.

Symmetry options

• :auto detects matrix symmetry. A symmetric matrix is tried as positive definite (SuiteSparse and Apple Accelerate only when its diagonal is positive) and refactored as symmetric indefinite when that fails. On SuiteSparse and Apple Accelerate, a singular matrix with a zero pivot throws a SingularException; a nearly singular one may return a large solution without an error. So may a singular one on Apple Accelerate, whose ordering can round the zero pivot of a pure-Neumann Laplacian to about eps. See the docstring of each backend's *_factorize.

• :spd, :definite, or 1 selects symmetric positive definite.

• :symmetric or 2 selects general symmetric.

• :unsymmetric or 0 selects general unsymmetric.

Ignored by :sparspak, which always factors as general unsymmetric LU. It too throws a SingularException for a singular matrix, with index 0 rather than a pivot position.

kwargs go to the chosen backend's *_factorize, which throws on a keyword it does not take. :sparspak takes none and throws an ArgumentError.

A SparseMatrixCSR (SparseMatricesCSR.jl) is also accepted: converted to SparseMatrixCSC first (no backend has a native CSR solve path – see docs/src/internals/csr_solvers.md), then factored exactly as above.

The form method refuses symmetrize with an ArgumentError: it also changes the load vector. Call A, F = assemble(a, l; symmetrize = true) and factor the matrix A.

Examples

using Bramble: sparse_factorize
fact = sparse_factorize(A; solver = :suitesparse, sym = :spd)
u = fact \ F

See also refactor!, pde_solve, suitesparse_factorize, suitesparse_qr_factorize, accelerate_factorize, mumps_factorize, sparspak_factorize.

refactor!(fact::Factorization, A::SparseMatrixCSC) -> Factorization
refactor!(fact::Factorization, a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing) -> Factorization

Recompute the numerical values of fact for updated matrix A (or assembled bilinear form a) reusing the existing symbolic factorization (fill-reducing ordering and elimination tree). A must have the exact same sparsity pattern as the matrix originally factored.

Dispatches automatically via multiple dispatch to the appropriate backend.

• SuiteSparseFactorization updates CHOLMOD or UMFPACK numeric values.

• AccelerateFactorization updates Apple Accelerate libSparse numeric values.

• MUMPSFactorization updates MUMPS multifrontal numerical factorization (job = 2).

• SparspakFactorization updates Sparspak's numeric LU values.

See also sparse_factorize, pde_solve.

SuiteSparseFactorization{T} <: Factorization{T}

Wrapper type representing a factorized SuiteSparse linear system (CHOLMOD Cholesky or UMFPACK LU), supporting in-place solves via LinearAlgebra.ldiv!, back-substitution via \, and non-allocating numeric refactoring via suitesparse_refactor!.

suitesparse_factorize(A::SparseMatrixCSC; sym = :auto, kwargs...) -> SuiteSparseFactorization
suitesparse_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing,
                      sym = :auto, kwargs...) -> SuiteSparseFactorization

Compute the sparse direct factorization of A (or the assembled matrix of a) using SuiteSparse (CHOLMOD for symmetric positive-definite systems, UMFPACK for unsymmetric).

Symmetry options

• :auto (default) detects symmetry. A symmetric A with a positive diagonal is tried with CHOLMOD sparse Cholesky, and when that fails (A is not positive definite) UMFPACK sparse LU factors it instead, so a symmetric indefinite matrix is solved and a singular one throws a SingularException. Any other A uses UMFPACK sparse LU. Under :auto the kwargs reach whichever factorisation runs, so a keyword only Cholesky takes (such as shift) errors after the fallback.

• :spd, :definite, or 1 selects symmetric positive definite (CHOLMOD Cholesky).

• :symmetric or 2 selects symmetric factorization.

• :unsymmetric or 0 selects general unsymmetric (UMFPACK LU).

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

Requires SuiteSparse.jl; call using SuiteSparse before calling this function.

Examples

using Bramble: suitesparse_factorize
using SuiteSparse

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0, symmetrize = true)
fact = suitesparse_factorize(A; sym = :spd)
u = fact \ F

See also suitesparse_solve, suitesparse_refactor!, pde_solve, assemble.

suitesparse_solve(A::SparseMatrixCSC, F::AbstractVector; sym = :auto, kwargs...) -> Vector
suitesparse_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
                  symmetrize = false, sym = :auto, kwargs...) -> VectorElement

Directly solve A u = F (or assemble(a, l) system) using SuiteSparse factorization. sym is as in suitesparse_factorize: symmetrize = true does not choose Cholesky, :auto still lets CHOLMOD decide.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

Requires SuiteSparse.jl; call using SuiteSparse before calling this function.

See also suitesparse_factorize, refactor!, pde_solve.

suitesparse_refactor!(fact::SuiteSparseFactorization, A::SparseMatrixCSC) -> SuiteSparseFactorization

Recompute the numeric factorization of A inside fact reusing the existing symbolic factorization (fill-reducing ordering and elimination tree). A must have the exact same sparsity pattern as the matrix originally factored.

See also refactor!, suitesparse_factorize.

suitesparse_qr_factorize(A::SparseMatrixCSC; tol = ..., ordering = ..., kwargs...) -> QRSparse
suitesparse_qr_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing,
                         kwargs...) -> QRSparse

Compute the sparse direct QR factorization of A (or the assembled matrix of a) using SuiteSparse's SPQR, suited to overdetermined least-squares systems and the rectangular blocks of a constrained saddle-point form. Unlike suitesparse_factorize, A need not be square, and this needs only SparseArrays – already a dependency of Bramble, no using SuiteSparse required.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

kwargs (tol, ordering, ...) are forwarded to SparseArrays.SPQR.qr.

Examples

using Bramble: suitesparse_qr_factorize
A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
fact = suitesparse_qr_factorize(A)
u = fact \ F

See also suitesparse_qr_solve, suitesparse_factorize, pde_solve.

suitesparse_qr_solve(A::SparseMatrixCSC, F::AbstractVector; tol = ..., ordering = ..., kwargs...) -> Vector
suitesparse_qr_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
                    kwargs...) -> VectorElement

Directly solve A u = F (exactly if square, least-squares if overdetermined) using SuiteSparse's SPQR sparse QR factorization. Needs only SparseArrays; no using SuiteSparse required. A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

See also suitesparse_qr_factorize, pde_solve.

AccelerateFactorization{T} <: Factorization{T}

Wrapper type representing a factorized Apple Accelerate (libSparse) linear system, supporting in-place solves via LinearAlgebra.ldiv!, back-substitution via \, and non-allocating symbolic reuse via accelerate_refactor!. Both solves take a vector or a matrix right-hand side of any real or complex element type, converted to the factorization's; an in-place ldiv! needs a destination that can hold the solution. Available only on macOS.

accelerate_factorize(A::SparseMatrixCSC; sym = :auto, kind = :auto) -> AccelerateFactorization
accelerate_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing,
                     sym = :auto, kind = :auto) -> AccelerateFactorization

Compute the sparse direct factorization of A (or the assembled matrix of a) using Apple Accelerate's native libSparse on macOS.

Symmetry and factorization options

• sym = :auto (default) detects symmetry. A Hermitian A (for a real A, a symmetric one) with a positive real diagonal is tried with Cholesky (SparseFactorizationCholesky), and when that fails (A is not positive definite) it is refactored as \(LDL^T\) (SparseFactorizationLDLT) on the same object. A Hermitian A with a non-positive diagonal uses \(LDL^T\) directly. Any other A, a complex symmetric one that is not Hermitian included, uses threshold partial pivoting LU (SparseFactorizationLUTPP). None of these reports a singular matrix itself, so :auto (and nothing) solves once against a fixed vector and throws a SingularException when the backward error is large. A singular A with a zero pivot throws, whatever the right-hand side. A nearly singular one may return a large solution without an error, and so may a singular one whose zero pivot rounds to about eps, as a pure-Neumann Laplacian's does. A structurally singular A (a zero row or column, say), which libSparse refuses itself, throws a SingularException too, for every sym and kind. The index of an Accelerate SingularException is always 0 and is not a pivot position: libSparse reports none, so Julia's "zero pivot found at index 0" message does not mean the first pivot.

• sym = :spd, :definite, or 1 (or kind = :cholesky) selects Hermitian positive definite Cholesky.

• sym = :symmetric or 2 (or kind = :ldlt) selects Hermitian indefinite \(LDL^T\). These two take a Hermitian A only: any other A, real or complex, throws an ArgumentError.

• sym = :unsymmetric or 0 (or kind = :lu / :lutpp) selects general unsymmetric LU with threshold partial pivoting.

• kind = :qr selects sparse QR factorization.

On a matrix, any keyword other than sym and kind throws a MethodError. The form method also takes dirichlet and dirichlet_components, and refuses symmetrize with an ArgumentError.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first, so a dense-backed form uses libSparse; accelerate_factorize(::Matrix) itself stays dense LAPACK.

Requires macOS and AppleAccelerate.jl; call using AppleAccelerate before calling this function.

Examples

using Bramble: accelerate_factorize
using AppleAccelerate

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0, symmetrize = true)
fact = accelerate_factorize(A; sym = :spd)
u = fact \ F

See also accelerate_solve, accelerate_refactor!, pde_solve, assemble.

accelerate_factorize(A::AbstractMatrix; sym = :auto, kind = :auto) -> LinearAlgebra.Factorization

Dense factorization of A, via ordinary LinearAlgebra.lu, LinearAlgebra.cholesky, or LinearAlgebra.qr – not Apple Accelerate's libSparse used by the SparseMatrixCSC method above, and not gated on macOS.

AppleAccelerate.jl (checked against v0.7.0) defines no dense lu, cholesky, getrf, potrf, gemm, or other LAPACK/BLAS bindings for Bramble to call. What its __init__ does is register Accelerate with libblastrampoline (BLAS.lbt_forward), so plain LinearAlgebra.lu/cholesky/qr are already running on Accelerate's BLAS/LAPACK the moment using AppleAccelerate has been evaluated, on every call site in Bramble or anywhere else, with no Bramble dispatch involved. This method adds no acceleration beyond that; it exists so dense callers can spell accelerate_factorize with the same sym/kind vocabulary as the sparse methods above.

Symmetry and factorization options

• sym = :auto (default) uses LinearAlgebra.cholesky if A is symmetric positive definite, otherwise LinearAlgebra.lu.

• sym = :spd, :definite, or 1 (or kind = :cholesky) uses LinearAlgebra.cholesky.

• sym = :unsymmetric or 0 (or kind = :lu) uses LinearAlgebra.lu.

• kind = :qr uses LinearAlgebra.qr.

• sym = :symmetric, 2, or kind = :ldlt is not implemented here. The dense analogue of a sparse LDLᵀ factorization is LinearAlgebra.bunchkaufman, which this method does not wrap; call it directly.

See also accelerate_solve.

accelerate_solve(A::SparseMatrixCSC, F::AbstractVector; sym = :auto, kind = :auto) -> Vector
accelerate_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
                 symmetrize = false, sym = :auto, kind = :auto) -> VectorElement

Directly solve A u = F (or assemble(a, l) system) using Apple Accelerate direct factorization. sym is as in accelerate_factorize: symmetrize = true does not choose Cholesky, :auto still tries it and falls back to \(LDL^T\). Any keyword other than sym and kind throws a MethodError.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first, so a dense-backed form uses libSparse; accelerate_factorize(::Matrix) itself stays dense LAPACK.

Requires macOS and AppleAccelerate.jl; call using AppleAccelerate before calling this function.

See also accelerate_factorize, refactor!, pde_solve.

accelerate_refactor!(fact::AccelerateFactorization, A::SparseMatrixCSC) -> AccelerateFactorization

Recompute the numeric factorization of A stored in fact reusing the existing symbolic factorization (fill-reducing ordering and sparsity analysis). A must have the exact same sparsity pattern as the matrix originally factored.

A failed refactor of any kind leaves the old factors unusable. Every later solve then throws the exception the refactor threw, until a refactor succeeds: a SingularException (its index is not a pivot position, as for accelerate_factorize) for a singular A with a zero pivot under sym = :auto or a structurally singular one, libSparse's error for an A that is not positive definite under an explicit sym = :spd. When :auto chose Cholesky and A is not positive definite, A is factored afresh as \(LDL^T\), as accelerate_factorize would, without the symbolic reuse. A Cholesky or \(LDL^T\) factorization, from :auto or chosen, refuses an A that is not Hermitian with an ArgumentError and keeps its old factors, as every kind does for an A with another number of stored entries.

See also refactor!, accelerate_factorize.

MUMPSFactorization{T} <: Factorization{T}

Wrapper type representing a factorized MUMPS linear system, supporting in-place solves via LinearAlgebra.ldiv!, back-substitution via \, and non-allocating reuse for transient PDE time stepping.

mumps_factorize(A::SparseMatrixCSC; sym = :auto, icntl = nothing, cntl = nothing) -> MUMPSFactorization
mumps_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing,
                sym = :auto, icntl = nothing, cntl = nothing) -> MUMPSFactorization

Compute the sparse direct multifrontal factorization of A (or the assembled matrix of a) using MUMPS.jl.

Symmetry options

• :auto: the default, which detects symmetry. A symmetric A is factored as symmetric positive definite (sym = 1, no pivoting), and when MUMPS reports a zero or negative pivot (INFOG(1) = -10 or INFOG(12) > 0) it is factored again as general symmetric \(LDL^T\) (sym = 2). Any other A uses unsymmetric LU (sym = 0).

• :spd, :definite, 1: symmetric positive definite (Cholesky / \(LL^T\)).

• :symmetric, 2: general symmetric (\(LDL^T\) with Bunch-Kaufman pivoting).

• :unsymmetric, 0: general unsymmetric LU.

A failed factorisation (INFOG(1) < 0) throws for every sym: a SingularException for a zero pivot (INFOG(1) = -10), an error otherwise.

The factorization's MUMPS memory is freed once it is unreachable; finalize(fact) frees it at once, after which fact must not be used.

Control parameters

• icntl: Optional dictionary or collection of pairs of integer control parameters (e.g., 7 => 1 for user/METIS ordering, 14 => 30 for memory relaxation).

• cntl: Optional dictionary or collection of pairs of real control parameters (e.g., 1 => 0.01 for numerical pivoting threshold).

On a matrix, any other keyword throws a MethodError. The form method also takes dirichlet and dirichlet_components, and refuses symmetrize with an ArgumentError.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

Requires MUMPS.jl; call using MUMPS before calling this function.

Examples

using Bramble: mumps_factorize
using MUMPS

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0, symmetrize = true)
fact = mumps_factorize(A; sym = :spd)
u = fact \ F

See also pde_solve, mumps_solve, assemble.

mumps_solve(A::SparseMatrixCSC, F::AbstractVector; sym = :auto, icntl = nothing, cntl = nothing) -> Vector
mumps_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
            symmetrize = false, sym = :auto, icntl = nothing, cntl = nothing) -> VectorElement

Directly solve A u = F (or assemble(a, l) system) using MUMPS direct factorization.

Symmetry options

• :auto: the default, as in mumps_factorize: symmetrize = true does not choose positive definite, and a symmetric indefinite A is refactored as \(LDL^T\).

• :spd, :definite, 1: symmetric positive definite.

• :symmetric, 2: general symmetric.

• :unsymmetric, 0: general unsymmetric.

icntl and cntl are as in mumps_factorize. Any other keyword throws a MethodError.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

Requires MUMPS.jl; call using MUMPS before calling this function.

Examples

using Bramble: mumps_solve
using MUMPS

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
u = mumps_solve(A, F)

See also mumps_factorize, refactor!, pde_solve.

mumps_refactor!(fact::MUMPSFactorization, A::SparseMatrixCSC) -> MUMPSFactorization

Recompute the numeric factorization of A inside fact reusing the existing symbolic factorization (fill-reducing analysis and ordering). A must have the exact same sparsity pattern as the matrix originally factored.

See also refactor!, mumps_factorize.

SparspakFactorization{T} <: Factorization{T}

Wrapper type representing a factorized Sparspak.jl sparse LU system, supporting in-place solves via LinearAlgebra.ldiv!, back-substitution via \, and non-allocating reuse for transient PDE time stepping.

sparspak_factorize(A::SparseMatrixCSC) -> SparspakFactorization
sparspak_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing) -> SparspakFactorization

Compute the pure-Julia sparse direct LU factorization of A (or the assembled matrix of a) using Sparspak.jl.

Sparspak has no binary dependency, so it factors matrices whose entries are not Float64/ComplexF64 – Float32, BigFloat, or a ForwardDiff.Dual – where SuiteSparse and MUMPS cannot. pde_solve does not hand it a Dual matrix: it factors the Float64 values of that matrix once and solves each partial against them (see pde_solve), so a Dual factor is built only by a direct call.

A singular matrix throws a SingularException, as it does for SuiteSparse, MUMPS and Accelerate, from the factorization or from the solve: a rank-deficient or structurally singular matrix at factorization, and a stored zero pivot or an empty row or column when a solve returns a non-finite entry for a finite right-hand side. Sparspak's own failure is a misleading ErrorException or an Inf in the solution. The index of the SingularException is 0, not a pivot position.

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

Requires Sparspak.jl; call using Sparspak before calling this function.

Examples

using Bramble: sparspak_factorize
using Sparspak

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
fact = sparspak_factorize(A)
u = fact \ F

See also sparspak_solve, sparspak_refactor!, pde_solve, assemble.

sparspak_solve(A::SparseMatrixCSC, F::AbstractVector) -> Vector
sparspak_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
               symmetrize = false) -> VectorElement

Directly solve A u = F (or assemble(a, l) system) using Sparspak's pure-Julia sparse direct LU factorization.

A singular matrix throws a SingularException with index 0 (see sparspak_factorize).

A BilinearForm on a dense or CSR backend is converted to SparseMatrixCSC first.

Requires Sparspak.jl; call using Sparspak before calling this function.

Examples

using Bramble: sparspak_solve
using Sparspak

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
u = sparspak_solve(A, F)

See also sparspak_factorize, sparspak_refactor!, pde_solve.

sparspak_refactor!(fact::SparspakFactorization, A::SparseMatrixCSC) -> SparspakFactorization

Recompute the numeric factorization of A inside fact reusing the existing symbolic factorization (fill-reducing ordering). A must have the exact same sparsity pattern as the matrix originally factored. A singular A throws a SingularException with index 0.

See also refactor!, sparspak_factorize.