Sparse linear algebra

SparseMatrixCSC(I, J, V, m, n)

Construct an m×n sparse matrix from triplet lists I (rows), J (columns) and V (values). Duplicate entries are summed.

SparseMatrixCSC(A::AbstractMatrix)

Convert a dense (or otherwise indexable) matrix to a SparseMatrixCSC, dropping exact zeros.

SparseMatrixCSR(I, J, V, m, n)

Construct an m×n CSR matrix from triplet lists I (rows), J (columns) and V (values). Duplicate entries are summed.

SparseMatrixCSR(A::AbstractMatrix)

Convert a dense (or otherwise indexable) matrix to a SparseMatrixCSR, dropping exact zeros.

analyze(A; ordering = :amd, perm = nothing) -> CholeskyAnalysis

Perform the symbolic phase (ordering, permutation, elimination tree, symbolic factorization, supernode detection) once, for reuse across multiple numeric factorizations with the same sparsity pattern.

amd(A; dense = 10.0, aggressive = true, variant = :faithful) -> perm

Compute a fill-reducing ordering of the symmetric structure of A.

Keyword arguments

• dense: dense-row threshold. Rows with more than max(dense·√n, 16) off-diagonal entries are removed and placed last (set dense < 0 to only remove completely dense rows).
• aggressive: enable aggressive absorption (element absorption and supervariable amalgamation) in the :faithful variant.
• variant: :faithful (port of the SuiteSparse AMD algorithm, the default), :classic (native quotient-graph minimum degree), or :hybrid (run both, plus the non-aggressive faithful variant, and return the ordering with the smallest estimated fill — ~3× more expensive than a single variant).

nested_dissection(A; min_size = 32) -> perm

Nested-dissection fill-reducing ordering for the symmetric structure of A (via the recursive level-structure separator of George & Liu). Places a vertex separator last at each recursion level.

Refactorization(A; method = :cholesky, ...) -> R

Cache the symbolic analysis (ordering + symbolic factor pattern) of the sparse matrix A once, so a sequence of matrices sharing A's sparsity pattern — the Newton steps of an optimization loop, the stiffness matrices across topology- optimization iterations, a time-stepping sequence — can be re-factorized paying only the numeric phase.

R = Refactorization(A; method = :lu, ordering = :amd)
u = R(A_new) \ b          # factorize A_new reusing the ordering, then solve
F = R(A_new)               # just the numeric factorization

method:

• :cholesky / :direct (default) — SPD; analyze + factorize.
• :ldlt — symmetric indefinite; analyze + ldlt.
• :lu — unsymmetric; lu_analyze + lu_factorize.
• :klu — unsymmetric block-triangular form; klu_analyze + klu_factorize.

Keyword arguments passed on construction (ordering, perm, supernode_min, relax, relax_fill) are forwarded to the symbolic phase; keyword arguments passed on each call (check, tol, shift, parallel, scale, method, …) are forwarded to the numeric phase. A_new must have the same sparsity pattern as A (the standard analyze/factorize contract); only the values may differ.

perm(F) -> Vector{Int}

The permutation associated with a factorization F (identity for a SparseTriangular factor, which carries no fill-reducing permutation).

factor_nnz(F) -> Int

The number of stored nonzeros in a factorization F (e.g. the sum of its factors).

cg(A, b; M=nothing, atol=0.0, rtol=1e-8, maxiter=nothing) -> (; x, iters, history)

Preconditioned conjugate gradient (Hestenes–Stiefel) for SPD A.

gmres(A, b; restart=min(20, n), M=nothing, Pr=nothing, atol=0.0, rtol=1e-8, maxiter=nothing)

Restarted GMRES with modified Gram–Schmidt Arnoldi, least-squares solved via a QR of the Hessenberg matrix. M is the left preconditioner (M⁻¹ A x = M⁻¹ b) and Pr an optional right preconditioner (two-sided M⁻¹ A Pr⁻¹ (Pr x) = M⁻¹ b).

bicgstab(A, b; M=nothing, atol=0.0, rtol=1e-8, maxiter=nothing)

BiCGStab (van der Vorst) for nonsymmetric A, with a left preconditioner M.

jacobi(A; ω = 1.0, shift = 0.0) -> Jacobi

Construct the (damped) Jacobi preconditioner M = D/ω from the diagonal D of the square matrix A. Applying it computes M \ b = ω D⁻¹ b elementwise. shift = ε adds ε to the diagonal first (D → D + εI), so the preconditioner stays well-defined even for a matrix with zero diagonal entries.

ssor(A; ω = 1.0, shift = 0.0) -> SSOR

Construct the SSOR preconditioner M = 1/(ω(2-ω)) · (D + ωL) D⁻¹ (D + ωU) of the square matrix A (D diagonal, L strict lower, U strict upper). Applying it computes M \ b with two triangular solves. ω ∈ (0, 2) is the relaxation factor (default 1.0 = symmetric Gauss–Seidel). shift = ε adds ε to the diagonal first (D → D + εI), so the preconditioner stays well-defined even for a matrix with zero diagonal entries.

ic(A; method = :ic0, τ = 0.0, k = 0, ordering = :amd, perm = nothing, check = true, shift = 0.0) -> preconditioner

Incomplete Cholesky factorization of the symmetric positive-definite matrix A (A ≈ L Lᴴ), applied with M \ b or ldiv!(M, b).

method

• :ic0 — zero fill-in (IC(0)): L shares A's lower pattern.
• :ict — threshold with the full symbolic fill pattern (ICT): AMD ordering + the complete Cholesky fill pattern, dropping entries below τ.
• :ick — level-k fill (IC(k)): keep entries whose shortest-path distance (with intermediate vertices < min(i,j)) is ≤ k + 1.

τ

Magnitude drop tolerance (off-diagonal entries of L with magnitude below τ are dropped; τ = 0 keeps the full pattern).

ordering (:ict only)

• :amd (default) — approximate minimum degree.
• :natural — natural order.
• perm — an explicit symmetric permutation (overrides ordering).

check

Mirrors cholesky: true throws PosDefException on a non-positive pivot.

shift

• ε ≥ 0 — factor A + εI instead of A (diagonal regularization, so M approximates (A + εI)⁻¹); 0 (default) means no shift.

ilu(A; method = :ilu0, τ = 0.0, k = 0, ordering = :amd, perm = nothing, shift = 0.0) -> preconditioner

Incomplete LU factorization of the square matrix A (A ≈ L U), applied with M \ b or ldiv!(M, b).

method

• :ilu0 — zero fill-in (ILU(0)): L/U share A's pattern, no pivoting.
• :ilut — threshold with fill (ILUTP): a fill-reducing column ordering (ordering) is applied, then Gilbert–Peierls partial pivoting with magnitude dropping (τ).
• :iluk — level-k fill (ILU(k)): keep entries whose shortest directed-graph distance is ≤ k + 1, no pivoting.

τ

Magnitude drop tolerance (off-diagonal entries below τ times the pivot are dropped). τ = 0 keeps the full pattern.

ordering (:ilut only)

• :amd (default) — column approximate minimum degree (colamd).
• :natural — natural column order.
• perm — an explicit column permutation (overrides ordering).

shift

• ε ≥ 0 — factor A + εI instead of A (diagonal regularization, so M approximates (A + εI)⁻¹); 0 (default) means no shift.

amg(A; method = :classical, strength = :classical, coarsening = :rs,
    interpolation = :modified, aggregation = :greedy,
    prolongation = :smoothed, smoother = :gs, cycle = :v, presmooth = 1,
    postsmooth = 1, θ = 0.25, max_levels = 20, coarse_size = 16,
    coarse_solver = :auto, ω = 2/3, max_agg_size = 8, near_nullspace = nothing,
    blocksize = 1, non_galerkin = false, ng_threshold = 0.0,
    strength_iters = 4, candidates = 2, air_degree = 1,
    aggressive = false, aggressive_levels = 1, trunc_factor = 0.0,
    max_elements = 0, smoother_tol = 0.0, smoother_level = 0,
    coordinates = nothing, shift = 0.0)

Build an algebraic multigrid preconditioner for the square matrix A.

Families (method)

• :classical — Ruge–Stüben AMG (C/F coarsening + interpolation).
• :smoothed_aggregation / :sa — smoothed aggregation.
• :rootnode — root-node aggregation AMG (MIS aggregation + energy-minimizing prolongation, near-nullspace aware).
• :air — approximate ideal restriction AMG (nonsymmetric; one-point interpolation + local ideal restriction).
• :mgr — multigrid reduction AMG.
• :adaptive / :alpha_sa — adaptive smoothed aggregation.
• :bootstrap — bootstrap AMG (adaptive near-nullspace via the hierarchy).

Components

• strength: :auto, :classical, :symmetric, :evolution, :algebraic_distance, :affinity. :auto picks :symmetric for the aggregation families (SA/root-node/adaptive/bootstrap, matching PyAMG/MueLu) and :classical for the C/F families (classical/AIR/MGR, matching hypre).
• coarsening (classical): :rs, :pmis, :hmis, :cljp, :falgout, :cr, :cgc, :cgce.
• interpolation (classical): :direct, :classical (unmodified Ruge–Stüben), :modified (default; matches PyAMG/hypre "modified classical"), :extended, :extended_i.
• aggregation (SA/root-node): :greedy/:coupled, :mis, :pairwise, :uncoupled (block-respecting), :uncoupled_mis, :structured (coordinate bricks), :zoltan (recursive coordinate bisection).
• prolongation (SA/root-node): :tentative, :smoothed, :energy_minimizing.
• smoother: :jacobi, :gs/:symmetric_gs, :forward_gs, :backward_gs, :l1_gs, :richardson, :chebyshev, :ilu, :ilut, :iluk, :block_jacobi, :block_gs/:block_sgs, :block_forward_gs, :block_backward_gs.
• cycle: :v, :w, :f.

Other options

• near_nullspace — an n×k matrix (or length-n vector) of near-nullspace vectors consumed by SA/root-node; defaults to the block constant/unit modes implied by blocksize.
• blocksize — the number of consecutive unknowns per node (block AMG and the block/uncoupled aggregation and block smoothers).
• non_galerkin / ng_threshold — drop small coarse-operator entries.
• aggressive / aggressive_levels — apply hypre-style aggressive (multipass) coarsening — PMIS/HMIS over the strength-of-strength graph — to the first aggressive_levels levels, automatically pairing it with extended+i (distance-2) interpolation.
• trunc_factor / max_elements — interpolation truncation (drop entries below trunc_factor·row-max) and per-row element cap (keep the largest max_elements), the hypre interp_trunc_factor / interp_max_elements knobs.
• smoother_tol / smoother_level — the drop tolerance / fill level for the :ilut / :iluk (Euclid-style ILU) smoothers.
• coordinates — an n×d matrix of node coordinates for :structured and :zoltan aggregation.
• shift — a nonnegative diagonal regularization ε: the hierarchy is built on A + εI (the shift propagates to every Galerkin coarse operator, so the coarse solves are well-posed even for semidefinite or mildly indefinite A). Defaults to 0.0 (no shift).

Defaults

The defaults follow the established parameterizations of PyAMG, hypre (BoomerAMG), and MueLu: θ = 0.25 (classical strength threshold), ω = 2/3 (damped Jacobi), a V(1,1) cycle with a symmetric Gauss–Seidel smoother, 20 max levels, a 16-unknown direct coarse solve, and 8 as the aggregation size cap. coarse_solver = :auto uses Bunch–Kaufman LDLᵀ for self-adjoint operators and LU otherwise.

AMGPreconditioner

An assembled algebraic multigrid preconditioner: an ordered list of levels (coarsest last), reusable scratch buffers, and (optionally) the rebuild settings for update!.

Fields

• levels: the Vector{AMGLevel{T}}, coarsest last.
• n: the matrix dimension.
• cycle: :v, :w, or :f.
• ws: per-level transient buffers (levels 1..L-1).
• config: rebuild settings for update!, or nothing if not reusable.

operator_complexity(M::AMGPreconditioner) -> number

Operator complexity: the total nonzeros across all levels relative to the fine operator (a standard AMG efficiency metric).