Linear Solvers

AbstractLinearSolver

Choice of linear solver for the state, adjoint and Dirichlet systems: DirectSolver (projected sparse Cholesky) or BlockCGSolver.

DirectSolver(; backend = :auto)

Projected sparse Cholesky factorisation (projected_cholesky); refactorised numerically when the conductivity changes. backend: :auto, :cholmod or :ldl.

BlockCGSolver(; preconditioner = :amg, rtol = 1e-10, maxiter = 0)

Projected block conjugate gradients (pbcg!) with an :amg, :jacobi or :none preconditioner, or a preconditioner choice that needs the discretization, such as DCTPreconditioner(disc) on uniform rectangle grids. maxiter = 0 uses the default of pbcg!. The previous solution is the initial guess of the next solve (warm start across conductivity updates).

X, stats = pbcg(A, B; nullspace = nothing, grounding = nothing, M = nothing, kwargs...)

Allocating convenience wrapper around BlockCGWorkspace and pbcg!. B may be a vector or an n × s matrix (on the CPU or the GPU).

pbcg!(X, ws, A, B; M = nothing, rtol = √eps, atol = 0, maxiter = 10n, recompute_every = 50,
      rank_tol = √eps)

Solve A X = Π B with the projected block conjugate gradient method (Krylov.block_cg! of the Krylov.jl fork) in the workspace ws (see BlockCGWorkspace), in place in X (which is used as initial guess), and return a BlockCGStats.

Π = I - V Vᵀ is the orthogonal projector onto V⊥ = range(A). Every residual and every preconditioned residual is projected, so the iteration never leaves V⊥ (round-off and preconditioners that do not preserve V⊥ cannot pollute the solution, and inconsistent data are handled by solving the projected, consistent problem). At the end the component in V is fixed by the grounding condition Wᵀx = 0 of the workspace.

M is a preconditioner applied by ModularEIT.apply_preconditioner! (JacobiPreconditioner, AMGPreconditioner, DCTPreconditioner, PolarPreconditioner) or nothing. Column j has converged when its residual is below atol + rtol ‖Π bⱼ‖; converged columns are removed from the block (deflation), dependent search directions are dropped (relative eigenvalue threshold rank_tol), and every recompute_every iterations the updated residual is replaced by the true residual Π(B - AX) to limit round-off drift.

BlockCGWorkspace(A, B; nullspace = nothing, grounding = nothing)

Preallocated storage for pbcg! with size(B, 2) right-hand sides: a Krylov.BlockCgWorkspace of the Krylov.jl fork. All n × s buffers have the array type of B, so passing device arrays gives a device workspace.

• nullspace: basis of V = ker(A) as an n × k matrix or a vector. Default: the constant vector (pure Neumann problem). It does not need to be orthonormal. Pass an n × 0 matrix for a positive definite A (e.g. a Dirichlet system).
• grounding: linear functionals W (n × k or vector) fixing the component in V; the returned solution satisfies Wᵀx = 0. Default: W = V, i.e. the solution orthogonal to the null space (minimum Euclidean norm). For EIT with the boundary sum fixed to zero, pass the indicator vector of the boundary degrees of freedom (see boundary_grounding). WᵀV must be invertible.

Result information of pbcg!.

• converged: all columns reached their tolerance
• iterations: number of block iterations
• residuals: final relative residuals ‖Π(bⱼ - A xⱼ)‖ / ‖Π bⱼ‖ per column
• compatibility_defect: ‖(I - Π) B‖ / ‖B‖, the part of the right-hand side outside range(A) = V⊥ that was removed (should be ≈ 0 for a consistent Neumann problem)

boundary_grounding(n, dofs; weights = nothing)

Grounding functional w with wᵢ = 1 (or weights) on the degrees of freedom dofs, so that wᵀx = Σ_{i ∈ dofs} xᵢ = 0. Pass it as grounding to pbcg.

projected_cholesky(A; nullspace = nothing, grounding = nothing, nrhs = 1, to_device = identity)

Sparse Cholesky solver for A X = B where A (a SparseMatrixCSC) is symmetric positive semidefinite with null space V and positive definite on V⊥, i.e. an SPD operator on ℝⁿ/V. Returns a ProjectedCholesky that solves with F \ B or ldiv!(X, F, B) for a vector or a block of right-hand sides.

• nullspace: basis of V (n × k matrix or vector; default: constants, the pure Neumann case).
• grounding: functionals W fixing the component in V, the solution satisfies Wᵀx = 0 (default W = V: solution orthogonal to V). For EIT, boundary_grounding(n, boundary_dofs) makes the boundary values sum to zero; see boundary_grounding.
• nrhs: number of right-hand sides to preallocate buffers for (other sizes reallocate once).
• to_device: converts matrices/vectors to a device, e.g. device_converter(ROCArray). Without a device factorisation for the matrix type, the factorisation runs on the host and solves copy the right-hand sides (works on every GPUArrays backend). With CUDA.jl and CUDSS.jl loaded, to_device = x -> x isa SparseMatrixCSC ? CuSparseMatrixCSR(x) : CuArray(x) factorises and solves on the GPU with cuDSS.
• backend: :cholmod (CPU, Float64, supernodal), :ldl (CPU, LDLFactorizations.jl, any floating-point type) or :auto (default: CHOLMOD for Float64 on the CPU, LDLᵀ for other types, the device backend when to_device is given). See also projected_ldl.

Right-hand sides are projected onto range(A) = V⊥ first, so inconsistent data are solved in the least-squares sense. Use refactor! after the matrix values change (same pattern).

projected_ldl(A; kwargs...)

projected_cholesky with the LDLᵀ backend of LDLFactorizations.jl (backend = :ldl): pure Julia, works for Float32, Float64 and other floating-point types, CPU only.

ProjectedCholesky

Factorisation object returned by projected_cholesky. Solve with F \ B or ldiv!(X, F, B); update the numerical values with refactor!.

refactor!(F::ProjectedCholesky, A)

Numerical refactorisation for new values of A with the same sparsity pattern (e.g. a new conductivity σ in the same mesh). Reuses the symbolic analysis (ordering, elimination tree).

ldiv!(X, F::ProjectedCholesky, B)

Solve A X = Π B and ground the solution (Wᵀ X = 0). X and B are vectors or n × s matrices (on the device of F).

X, stats = pbminres(A, B; nullspace = nothing, grounding = nothing, scaling = :jacobi, kwargs...)

Allocating convenience wrapper around ProjectedMinresWorkspace and pbminres!.

stats = pbminres!(X, ws, A, B; atol = 0, rtol = √eps, itmax = 0)

Solve A X = Π B with Krylov.jl's block MINRES in the workspace ws (see ProjectedMinresWorkspace) and ground the result (Wᵀ X = 0). X, B are vectors or n × s matrices. Tolerances refer to the Frobenius norm of the (scaled) block residual, as in Krylov.jl. itmax = 0 uses Krylov.jl's default 2n/s.

Right-hand sides whose columns are linearly dependent (or zero) are handled by solving for an orthonormal basis of their span and recombining.

ProjectedMinresWorkspace(A, B; nullspace = nothing, grounding = nothing,
                         scaling = :jacobi, diagonal = nothing)

Preallocated storage for pbminres! with size(B, 2) right-hand sides (Krylov.jl block-MINRES workspace plus projection/grounding buffers, all created with similar(B, …)). nullspace and grounding work as in BlockCGWorkspace. scaling = :jacobi solves the symmetrically scaled system D^{-1/2} A D^{-1/2} with D = diag(A); pass diagonal if diag(A) is not available for the matrix type (e.g. on the GPU).

columnwise = true solves the right-hand sides one after another with Krylov.jl's single-vector MINRES instead of block MINRES. Use it on the GPU: Krylov.jl's block MINRES currently fails there because of a method ambiguity between cuBLAS.jl and GPUArrays.jl in triangular rdiv!.

Result information of pbminres!: converged, iterations, Krylov.jl status string and compatibility_defect = ‖(I - Π) B‖ / ‖B‖.

JacobiPreconditioner(A; to_device = identity)

Diagonal (Jacobi) preconditioner M = diag(A). A must be a host matrix; to_device converts the stored inverse diagonal (e.g. device_converter(CuArray)).

AMGPreconditioner(A, s; to_device = identity, sweeps = 1, max_coarse = 64, max_levels = 10)

Symmetric smoothed-aggregation AMG V-cycle for s right-hand sides at once.

The hierarchy (aggregation, prolongations P, Galerkin coarse operators PᵀAP) is built on the host with AlgebraicMultigrid.jl, using the constants as near-null space. The cycle itself only uses sparse × dense products, damped Jacobi smoothing and a dense pseudo-inverse on the coarsest level, so it runs unchanged on any GPU backend when to_device moves matrices there, e.g. to_device = device_converter(CuArray) (or ROCArray, oneArray, MtlArray); see device_converter.

With the same number of pre- and post-smoothing sweeps and R = Pᵀ the cycle is a symmetric operator, positive definite on the complement of the constants, as required by CG. The Jacobi damping is ω = 4 / (3 ρ(D⁻¹A)), with ρ estimated by power iteration.

apply_preconditioner!(Z, M, R)

Z ← M⁻¹ R for a block of vectors. M = nothing is the identity.

_gram!(G, X, Y)

G ← XᵀY for tall-skinny blocks (n × s with n ≫ s): Krylov.kgram! (one GEMV per column for double-precision blocks with 2–8 columns, where GEMM kernels are slow).

to_device = device_converter(ArrayType)

Conversion function for the to_device keyword of the solvers and preconditioners on any GPUArrays backend: sparse matrices become DeviceSparseMatrixCSR, dense arrays ArrayType(x). Examples:

device_converter(CuArray)      # NVIDIA (CUDA.jl)
device_converter(ROCArray)     # AMD (AMDGPU.jl)
device_converter(oneArray)     # Intel (oneAPI.jl)
device_converter(MtlArray)     # Apple (Metal.jl, Float32 only)
device_converter(JLArray)      # CPU reference implementation (JLArrays.jl), for testing

On NVIDIA GPUs, x -> x isa SparseMatrixCSC ? CuSparseMatrixCSR(x) : CuArray(x) uses cuSPARSE instead of the generic kernel and enables the cuDSS direct solver.

DeviceSparseMatrixCSR(A::SparseMatrixCSC, ArrayType)

Compressed sparse row matrix whose index and value arrays are stored with ArrayType (CuArray, ROCArray, oneArray, MtlArray, JLArray, or Array). Supports mul!(Y, A, X, α, β) for dense vectors/matrices on the same backend through a generic KernelAbstractions kernel. Usually created by device_converter.

DCTPreconditioner(disc; variant = :constant)

Choice of the DCT preconditioner for BlockCGSolver on a uniform rectangle grid of bilinear elements (see structured_grid): the exact inverse of the system for a constant conductivity, applied with fast cosine transforms (O(n log n) per application), for every electrode model and for current- and voltage-driven problems. The number of CG iterations does not grow with the mesh size, only with the conductivity contrast.

• variant = :constant: σ̄ K with the geometric mean σ̄ of the nodal conductivities; iterations grow like √(σmax/σmin).
• variant = :scaled: S^½ K S^½ with the nodal conductivities S (Concus–Golub): for smooth conductivities the iterations then hardly depend on the contrast. Not for jumps: there the iteration count grows with mesh refinement.

The setup costs one transform solve per contact node (complete electrode model) and per constrained node (voltage-driven problems), once per forward model.

dct_preconditioner(disc, fm; system = :neumann, variant = :constant)

DCT preconditioner for the current-driven system fm.A (system = :neumann) or the voltage-driven block fm.A_ff (:dirichlet) of a forward model, for use with pbcg! (M = …). Call update_preconditioner!(M, A) after every change of the matrix values.

update_preconditioner!(M, A)

Update a fast-transform preconditioner (from dct_preconditioner or polar_preconditioner) to new values of its system matrix A (same sparsity pattern).

StructuredGrid

A uniform tensor grid of nx × ny nodes with spacings hx, hy, and the permutation perm from lexicographic node order (x index fastest) to the u dofs of a discretization, with the data of the transform solver. Built by structured_grid(disc).

structured_grid(disc)

The StructuredGrid of a discretization on a uniform rectangle grid (bilinear u space), for the DCT preconditioner. Back end contract (optional).

dct_neumann_solve(grid, B)

K₂⁺ B for the Q1 Neumann stiffness matrix K₂ of the structured grid, with B (nx ny × s, or a vector) in lexicographic node order: the solution with zero trapezoidal mean (∫ u = 0); B must be orthogonal to the constants for K₂ X = B to hold.

AbstractFastPreconditioner

Preconditioner choice for BlockCGSolver that inverts the constant-conductivity system with fast transforms: DCTPreconditioner (uniform rectangle grids), PolarPreconditioner (rotationally symmetric disk meshes).

PolarPreconditioner(disc; variant = :constant)

Choice of the FFT preconditioner for BlockCGSolver on a rotationally symmetric disk mesh of linear triangles, e.g. polar_grid (Ferrite back end): the constant-conductivity system is inverted with an FFT in the angle and one tridiagonal solve per angular mode along the radius (O(n log n)), for every electrode model and for current- and voltage-driven problems. The radial spacing is arbitrary (e.g. graded towards the boundary). Variants and costs as for DCTPreconditioner.

On a conformally mapped mesh (conformal_grid (Ferrite back end)) pass the disk mesh as reference: the stiffness matrix of the disk then preconditions the one of the mapped domain; since the map is conformal, both are spectrally equivalent with constants close to 1.

polar_preconditioner(disc, fm; system = :neumann, variant = :constant)

FFT preconditioner (see PolarPreconditioner) for the current-driven system fm.A or the voltage-driven block fm.A_ff of a forward model on a disk mesh, for pbcg!; call update_preconditioner! after every change of the matrix values.

polar_grid(nr, nθ; radius = 1, boundary_spacing = nothing, offset = 0)

Triangle mesh of the regular nθ-gon inscribed in the circle of the given radius: a centre node and nr rings of nθ nodes at the angles offset + 2π j / nθ. The centre is connected to the first ring by a fan, and every ring-to-ring quadrilateral is split along the same diagonal, so the mesh is invariant under the rotation by 2π/nθ (as required by the FFT solver, see PolarPreconditioner). Rings are uniformly spaced, or geometrically graded towards the boundary with the outermost spacing boundary_spacing (EIT sensitivity and the electrode edges are concentrated at the boundary). For L electrodes from angular_electrodes with coverage c, choose nθ as a multiple of 2L/c (e.g. 4L for c = 1/2): then the electrode edges lie on nodes and all electrodes consist of the same number of boundary facets.

PolarStructure

Rotationally symmetric node structure of a disk mesh (nr rings of nθ nodes, optionally a centre node) with the permutation perm from structure order to the u dofs, and the mode factorisations of the FFT solver. Built by polar_structure(disc).

polar_structure(disc)

The PolarStructure of a discretization on a polar disk mesh, for the polar FFT preconditioner. Back end contract (optional).

fast_neumann_solve(structure, B)

Moore–Penrose pseudo-inverse of the stiffness matrix of a PolarStructure (or of a StructuredGrid) applied to B (vector or n × s, in structure order): FFT in the angle, one tridiagonal solve per angular mode.

ConformalMap(boundary; method = :auto, center = nothing, modes = nothing, max_modes = 4096,
             boundary_tol = 1e-8, tol = 1e-13, maxiter = 2000, relaxation = 1,
             boundary_derivative = nothing)

Conformal map Φ from the unit disk onto the simply connected domain bounded by boundary, normalised by Φ(0) = center and Φ'(0) > 0. boundary is a closed curve t ↦ (x, y) on [0, 2π) or a vector of polygon vertices; center must lie inside (default: the centroid if it does, otherwise the point farthest from the boundary).

• method = :theodorsen: star-shaped domains that are not too far from a disk (|d log ρ/dϑ| < 1 in polar coordinates about the centre); fixed-point iteration.
• method = :wegmann: any smooth Jordan domain (e.g. non-star-shaped); started from Symm's integral equation and polished by Newton's method (quadratic convergence). Uses boundary_derivative t ↦ (x', y') if given, else central differences. Corners are not resolved well.
• method = :auto: Theodorsen if it converges, Wegmann otherwise.

modes Fourier modes resolve the boundary correspondence. By default (modes = nothing) they are doubled from 256 up to max_modes until boundary_error ≤ boundary_tol. Domains where |Φ'| varies strongly along the boundary need many: this "crowding" grows exponentially with the aspect ratio (the tips of a 3:1 ellipse need ~4000 modes), so for elongated domains the disk is a poor model domain. An ArgumentError is thrown when no resolution succeeds; a warning when the tolerance is not reached.

Evaluate with Φ(w) (complex w, |w| ≤ 1) and map_derivative(Φ, w). Fields: center, coefficients (Taylor coefficients of log((Φ(w) - center)/w)), boundary_error (largest deviation of Φ(e^{iθ}) from the boundary between the nodes, relative to the size of the domain), iterations, method.

map_derivative(Φ::ConformalMap, w)

Derivative Φ'(w) of the conformal map; |Φ'| is the local length scale factor (boundary current densities and contact impedances transform with it).

ConformalGrid

A mesh grid of a domain Ω obtained by mapping the polar mesh reference of the unit disk with the conformal map map (same connectivity and numbering). See conformal_grid.

conformal_grid(Φ::ConformalMap, nr, nθ; boundary_spacing = nothing, offset = 0)

Mesh of the domain Φ(𝔻): the polar_grid of the unit disk with the nodes mapped by Φ. Since Φ is conformal, the triangles keep their shapes up to O(h) and the ring grading carries over to the boundary of Ω (scaled by |Φ'|). All computations (electrode models, phantoms, objectives) use the mapped mesh; the disk mesh reference provides the fast preconditioner, PolarPreconditioner(disc; reference = cg.reference).