diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl
index ca96764..12465fa 100644
--- a/src/MatrixCovers.jl
+++ b/src/MatrixCovers.jl
@@ -33,6 +33,7 @@ include("initializers.jl")   # the start menu; consumed by both solver families
 include("soft_covers.jl")
 include("powermean.jl")        # the PowerMean soft covers, the soft-cover default
 include("sparse_cholesky.jl")  # CHOLMOD factorizations for the minimal-cover solvers
+include("tree_cholesky.jl")    # fill-free factorization of forest-patterned matrices, for preconditioning
 include("minimal_covers.jl")
 include("active_set.jl")       # exact solution of the AbsLog{2} problem from a near-optimal point
 include("transversal.jl")      # transversal-tight minimal covers
diff --git a/src/active_set.jl b/src/active_set.jl
index e829ee9..e6fe98b 100644
--- a/src/active_set.jl
+++ b/src/active_set.jl
@@ -223,7 +223,7 @@ function _TreeLaplacian(F::_Forest{T}, supp; budgets::NamedTuple=(;)) where {T}
     else
         idx, L = _tree_laplacian_triplets(F, supp)
     end
-    fac = L isa SparseMatrixCSC{Float64,Int} ? _laplacian_factor(L; budgets...) : _laplacian_factor(L)
+    fac = L isa SparseMatrixCSC ? _laplacian_factor(L; budgets...) : _laplacian_factor(L)
     return _TreeLaplacian{T,typeof(fac)}(idx, fac, size(L, 1))
 end
 
@@ -346,20 +346,30 @@ function _tree_laplacian_triplets(F::_Forest{T}, supp) where {T}
     return idx, sparse(I, J, V, nfree, nfree)
 end
 
-# Jacobi-preconditioned CG on a pinned tree Laplacian whose sparse Cholesky
-# factor would exceed the fill or flop budget. Such Laplacians arise from
-# expander-like supports, which fill in under elimination but are well
-# conditioned, so CG converges in few iterations. A solve that CG does not
-# finish falls back to the factorization.
-struct _LaplacianCG
-    L::SparseMatrixCSC{Float64,Int}   # upper triangle, as analyzed
-    S::SparseMatrixCSC{Float64,Int}   # both triangles, for products
-    dg::Vector{Float64}
-    F::SparseCholesky
-    factored::Base.RefValue{Bool}
+# Jacobi-preconditioned CG, with iterative refinement, on a pinned tree
+# Laplacian: in `Float64`, one whose sparse Cholesky factor would exceed the
+# fill or flop budget; in other element types, every sparse one, since CHOLMOD
+# factors only `Float64`. Such Laplacians arise from expander-like supports,
+# which fill in under elimination but have clustered spectra, so CG converges
+# in few iterations. A solve that refinement does not finish falls back to a
+# factorization: the analyzed CHOLMOD factor `F` in `Float64`, a dense Cholesky
+# factorization otherwise.
+struct _LaplacianCG{T}
+    L::SparseMatrixCSC{T,Int}         # upper triangle
+    S::SparseMatrixCSC{T,Int}         # both triangles, for products
+    dg::Vector{T}
+    F::Union{SparseCholesky,Nothing}  # analyzed, `Float64` only
+    fallback::Base.RefValue{Any}      # the factorization once CG gave up; `nothing` before
     maxiter::Int                      # CG iterations per pass
+    err::Base.RefValue{T}             # estimated 2-norm forward error of the last solve
 end
 
+# Estimated 2-norm forward error of the last `_laplacian_solve!`. A
+# factorization is taken as exact at the certification floor of the finish.
+_solve_error(::SparseCholesky) = 0.0
+_solve_error(::Cholesky) = 0.0
+_solve_error(C::_LaplacianCG) = C.err[]
+
 function _laplacian_factor(L::SparseMatrixCSC{Float64,Int}; fillbudget::Real=LSQR_FILL_BUDGET,
                            flopbudget::Real=LSQR_FLOP_BUDGET, cgmaxiter::Int=10 * size(L, 1) + 100)
     F = SparseCholesky()
@@ -367,41 +377,91 @@ function _laplacian_factor(L::SparseMatrixCSC{Float64,Int}; fillbudget::Real=LSQ
     analyze!(F, L)
     if sizeof(Float64) * factor_entries(F) > fillbudget || factor_flops(F) > flopbudget * nnz(L)
         S = L + triu(L, 1)'
-        return _LaplacianCG(L, S, Vector(diag(S)), F, Ref(false), cgmaxiter)
+        return _LaplacianCG{Float64}(L, S, Vector(diag(S)), F, Ref{Any}(nothing), cgmaxiter, Ref(0.0))
     end
     factorize!(F, L)
     return F
 end
-# Wider types: dense Cholesky of the upper triangle.
-_laplacian_factor(L::SparseMatrixCSC) = _laplacian_factor(Matrix(L))
+# Element types CHOLMOD does not factor; the budgets do not apply.
+function _laplacian_factor(L::SparseMatrixCSC{T,Int}; fillbudget::Real=LSQR_FILL_BUDGET,
+                           flopbudget::Real=LSQR_FLOP_BUDGET, cgmaxiter::Int=10 * size(L, 1) + 100) where {T}
+    size(L, 1) == 0 && return _laplacian_factor(Matrix(L))
+    S = L + triu(L, 1)'
+    return _LaplacianCG{T}(L, S, Vector(diag(S)), nothing, Ref{Any}(nothing), cgmaxiter, Ref(zero(T)))
+end
 _laplacian_factor(L::Matrix) = LinearAlgebra.cholesky!(Symmetric(L, :U))
 
 _laplacian_solve!(x::Vector{Float64}, F::SparseCholesky, b::Vector{Float64}) =
     isempty(b) ? x : solve!(x, F, CHOLMOD_A, b)
 _laplacian_solve!(x, F::Cholesky, b) = ldiv!(x, F, b)
 
-function _laplacian_solve!(x::Vector{Float64}, C::_LaplacianCG, b::Vector{Float64})
-    if !C.factored[]
+function _laplacian_solve!(x::Vector{T}, C::_LaplacianCG{T}, b::Vector{T}) where {T}
+    if C.fallback[] === nothing
         n = length(b)
-        fill!(x, 0.0)
-        work = ntuple(_ -> zeros(n), 4)
-        # The multipliers are routed sums of the solution along the trees, so
-        # certification needs a normwise backward error of a few units of
-        # roundoff, as from a factorization. `2 max(dg)` bounds `‖S‖₂` for a
-        # Laplacian; the first pass estimates `‖x‖`.
+        fill!(x, zero(T))
+        rdd, r, z, d, Ad, e = ntuple(_ -> zeros(T, n), 6)
         Smul! = (y, v) -> mul!(y, C.S, v)
-        _pcg!(Smul!, x, C.dg, b, work..., C.maxiter, sqrt(eps(Float64)) * norm(b))
-        tol = 8 * eps(Float64) * (2 * maximum(C.dg) * norm(x) + norm(b))
-        _, ok = _pcg!(Smul!, x, C.dg, b, work..., C.maxiter, tol)
-        ok && return x
-        factorize!(C.F, C.L)
-        C.factored[] = true
+        η = sqrt(eps(T))
+        _, ok = _pcg!(Smul!, x, C.dg, b, r, z, d, Ad, C.maxiter, η * norm(b))
+        # Iterative refinement with the residual in compensated arithmetic.
+        # CG leaves a forward error of its normwise backward error times the
+        # condition number, and the multipliers, routed sums of the solution
+        # along the trees, need the solution to the certification floor. A
+        # correction solved to relative residual `η` shrinks the error by a
+        # factor of about `κ η`, so the size of the last correction estimates
+        # the error that remains. Refinement that stops contracting hands the
+        # system to the factorization.
+        δprev = T(Inf)
+        while ok
+            _residual_dd!(rdd, b, C.S, x)
+            fill!(e, zero(T))
+            _, ok = _pcg!(Smul!, e, C.dg, rdd, r, z, d, Ad, C.maxiter, η * norm(rdd))
+            ok || break
+            x .+= e
+            δ = norm(e)
+            if δ <= 8 * eps(T) * norm(x)
+                C.err[] = δ
+                return x
+            end
+            δ <= δprev / 100 || break
+            δprev = δ
+        end
+        C.fallback[] = C.F === nothing ? _laplacian_factor(Matrix(C.L)) : factorize!(C.F, C.L)
+    end
+    C.err[] = zero(T)
+    return _laplacian_solve!(x, C.fallback[], b)
+end
+
+# `r = b - S x` with every product and sum carried in twice the working
+# precision by error-free transformations, so that `r` is accurate to a few
+# units of roundoff in its own magnitude rather than in that of `b` and `S x`.
+# `fma` must be correctly rounded in `T`.
+function _residual_dd!(r::Vector{T}, b::Vector{T}, S::SparseMatrixCSC{T,Int}, x::Vector{T}) where {T}
+    axes(r) == axes(b) == axes(x) == (axes(S, 1),) == (axes(S, 2),) ||
+        throw(DimensionMismatch("residual, right-hand side, matrix, and solution must share axes"))
+    copyto!(r, b)
+    lo = zeros(T, length(r))
+    rv = rowvals(S)
+    nz = nonzeros(S)
+    for j in axes(S, 2)
+        xj = x[j]
+        for k in nzrange(S, j)
+            i = rv[k]
+            p = nz[k] * xj
+            pe = fma(nz[k], xj, -p)   # `p + pe` is the exact product
+            hi = r[i]
+            s = hi - p
+            t = s - hi
+            lo[i] += (hi - (s - t)) + (-p - t) - pe
+            r[i] = s
+        end
     end
-    return solve!(x, C.F, CHOLMOD_A, b)
+    r .+= lo
+    return r
 end
 
-# Per-tree values `y` minimizing `½‖L y - rhs‖` with pinned trees at zero; `rhs`
-# is indexed by tree.
+# Per-tree values `y` minimizing `½‖L y - rhs‖` with pinned trees at zero, and
+# the estimated 2-norm error of the solve; `rhs` is indexed by tree.
 function _tree_solve(K::_TreeLaplacian, rhs::Vector{T}) where {T}
     b = zeros(T, K.nfree)
     for (t, i) in enumerate(K.idx)
@@ -413,10 +473,11 @@ function _tree_solve(K::_TreeLaplacian, rhs::Vector{T}) where {T}
     for (t, i) in enumerate(K.idx)
         i > 0 && (y[t] = x[i])
     end
-    return y
+    return y, _solve_error(K.fac)
 end
 
-# Minimizer of the objective with every edge of `F` tight.
+# Minimizer of the objective with every edge of `F` tight, and the estimated
+# error of the Laplacian solve behind it.
 function _tight_minimizer(F::_Forest{T}, K::_TreeLaplacian, supp::EdgeList) where {T}
     rhs = zeros(T, F.ntree)
     for (e, (p, q)) in enumerate(supp.edges)
@@ -426,8 +487,8 @@ function _tight_minimizer(F::_Forest{T}, K::_TreeLaplacian, supp::EdgeList) wher
         rhs[tp] -= k
         rhs[tq] += k
     end
-    y = _tree_solve(K, rhs)
-    return [y[F.tree[v]] + F.off[v] for v in eachindex(F.tree)]
+    y, δ = _tree_solve(K, rhs)
+    return [y[F.tree[v]] + F.off[v] for v in eachindex(F.tree)], δ
 end
 
 function _tight_minimizer(F::_Forest{T}, K::_TreeLaplacian, s::_DiffGridQP) where {T}
@@ -438,8 +499,8 @@ function _tight_minimizer(F::_Forest{T}, K::_TreeLaplacian, s::_DiffGridQP) wher
     for v in eachindex(F.tree, g)
         rhs[F.tree[v]] -= g[v]
     end
-    y = _tree_solve(K, rhs)
-    return [y[F.tree[v]] + F.off[v] for v in eachindex(F.tree)]
+    y, δ = _tree_solve(K, rhs)
+    return [y[F.tree[v]] + F.off[v] for v in eachindex(F.tree)], δ
 end
 
 # Add `∑ₑ dₑ gₑ` to `g`, where `dₑ = v[p] - v[q] - c[e]` when `withc` and
@@ -544,16 +605,22 @@ vector of nonnegative weights indexed like the constraints (edge positions of
 an `EdgeList`, linear indices of a `DiffGrid`'s `C`), such as the multipliers of
 an augmented-Lagrangian iteration, the candidates are the constraints with
 positive weight or residual below `ptol`, largest weight first and then
-tightest first. `certified` reports that the returned `u` satisfies the KKT
-conditions to tolerances proportional to the largest magnitude in `c` and `u0`;
-when it is false, `u` is `u0`. `kkt = (; primal, ptol, dual, dtol)` gives the
-largest constraint violation of the returned `u` and the most negative active
-multiplier (as nonnegative magnitudes), with the tolerances they were certified
-against; `primal` and `dual` are `NaN` when `certified` is false. Each step
-solves a graph-Laplacian system, and `maxsteps` bounds their number. `budgets`, a `NamedTuple` with fields
-`fillbudget` and/or `flopbudget`, sets the limits above which these systems are
-solved by conjugate gradients rather than a sparse Cholesky factorization.
-The constant on each connected component of the graph is taken from `u0`.
+tightest first. Constraints on no cycle of the support graph are activated
+ahead of the candidates: they are tight with zero multiplier at the minimizer,
+so a solve error must not let them appear violated. `certified` reports that the
+returned `u` satisfies the KKT conditions to tolerances proportional to the
+largest magnitude in `c` and `u0`; when it is false, `u` is `u0`.
+`kkt = (; primal, ptol, dual, dtol)` gives the largest constraint violation of
+the returned `u` and the most negative active multiplier (as nonnegative
+magnitudes), with the tolerances they were certified against; `primal` and
+`dual` are `NaN` when `certified` is false. Each step solves a graph-Laplacian
+system, and `maxsteps` bounds their number. `budgets`, a `NamedTuple` with
+fields `fillbudget` and/or `flopbudget`, sets the limits above which these
+systems are solved by conjugate gradients rather than a sparse Cholesky
+factorization; the tolerances of the certificate then rise to twice the
+estimated error of the solve, so that `certified` never claims more accuracy
+than the solve delivers. The constant on each connected component of the graph
+is taken from `u0`.
 """
 function _polish_difference_qp(edges::Vector{Tuple{Int,Int}}, cvals::Vector{T}, u0::Vector{T}; kwargs...) where {T}
     return _polish_difference_qp(EdgeList{T}(edges, cvals, -oneunit(T)), u0; kwargs...)
@@ -581,6 +648,8 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{
                      _seed_candidates(u0, s, ptol))
         sort!(cand; by=e -> (-seed[e], _resid(u0, s, e)))
     end
+    # Constraints on no cycle come first; Kruskal then skips them if they recur.
+    cand = vcat(_acyclic_constraints(s, nV), cand)
     uf = collect(1:nV)
     find(x) = (while uf[x] != x; uf[x] = uf[uf[x]]; x = uf[x]; end; x)
     inW = _active_flags(s)
@@ -600,13 +669,13 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{
     r = zeros(T, E)
     nsteps = 0
     u = u0
-    u, nsteps = _dual_feasible!(inW, W, λ, s, nV, dtol, nsteps, maxsteps; budgets)
+    u, nsteps, ptol_eff, dtol_eff = _dual_feasible!(inW, W, λ, s, nV, ptol, nsteps, maxsteps; budgets)
     certified = false
     nokkt = (; primal=T(NaN), ptol, dual=T(NaN), dtol)
     while nsteps <= maxsteps
-        e, de = _most_violated(u, inW, s, -ptol)
+        e, de = _most_violated(u, inW, s, -ptol_eff)
         if e == 0
-            certified = all(k -> λ[k] >= -dtol, W)
+            certified = all(k -> λ[k] >= -dtol_eff, W)
             break
         end
         p, q = _edge(s, e)
@@ -635,7 +704,7 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{
             rhs = zeros(T, F.ntree)
             rhs[F.tree[p]] += one(T)
             rhs[F.tree[q]] -= one(T)
-            y = _tree_solve(K, rhs)
+            y, _ = _tree_solve(K, rhs)
             z = [y[F.tree[v]] for v in 1:nV]
             a = z[p] - z[q]
             a > 0 || return u0, false, nsteps, nokkt
@@ -663,14 +732,14 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{
         end
         added || break
         # Recompute from `W` alone, so rounding does not accumulate.
-        u, nsteps = _dual_feasible!(inW, W, λ, s, nV, dtol, nsteps, maxsteps; budgets)
+        u, nsteps, ptol_eff, dtol_eff = _dual_feasible!(inW, W, λ, s, nV, ptol, nsteps, maxsteps; budgets)
     end
     certified || return u0, false, nsteps, nokkt
     # Restore each component's constant from `u0`. (`u` is reassigned above, so
     # capturing it in a closure would box it; `ustar` is bound once.)
     ustar = u
-    kkt = (; primal=max(zero(T), -minimum(e -> _resid(ustar, s, e), 1:_nconstraints(s))), ptol,
-           dual=max(zero(T), -minimum(k -> λ[k], W; init=zero(T))), dtol)
+    kkt = (; primal=max(zero(T), -minimum(e -> _resid(ustar, s, e), 1:_nconstraints(s))), ptol=ptol_eff,
+           dual=max(zero(T), -minimum(k -> λ[k], W; init=zero(T))), dtol=dtol_eff)
     uf .= 1:nV
     _foreach_constraint(s) do p, q
         rp, rq = find(p), find(q)
@@ -693,20 +762,28 @@ function _deactivate!(inW, W::Vector{Int}, λ, w::Int)
 end
 
 # Drop the edges of `W` with negative multipliers until none remain, and return
-# the tight minimizer of `W` with the step count; `λ` holds the multipliers.
-function _dual_feasible!(inW, W::Vector{Int}, λ::Vector{T}, supp, nV::Int, dtol,
+# the tight minimizer of `W`, the step count, and the primal and dual
+# tolerances that minimizer can be certified against; `λ` holds the
+# multipliers. The tolerances start from `ptol` and `ptol` times the number of
+# constraints and rise with the error of the Laplacian solve: an error `δ` in
+# the per-tree values moves a residual by at most `2δ`, and a routed
+# multiplier, a sum of residuals over a subtree, by at most `2δ` per constraint.
+function _dual_feasible!(inW, W::Vector{Int}, λ::Vector{T}, supp, nV::Int, ptol,
                          nsteps::Int, maxsteps::Int; budgets::NamedTuple=(;)) where {T}
+    E = max(1, _nsupported(supp))
     while true
         F = _Forest(nV, supp, W)
         K = _TreeLaplacian(F, supp; budgets)
-        u = _tight_minimizer(F, K, supp)
+        u, δ = _tight_minimizer(F, K, supp)
+        ptol_eff = max(ptol, T(2δ))
+        dtol_eff = ptol_eff * E
         _route!(λ, _difference_gradient(u, supp), F, supp)
         nW = length(W)
         for e in W
-            λ[e] < -dtol && (inW[e] = false)
+            λ[e] < -dtol_eff && (inW[e] = false)
         end
         filter!(e -> inW[e], W)
-        (length(W) < nW && nsteps < maxsteps) || return u, nsteps
+        (length(W) < nW && nsteps < maxsteps) || return u, nsteps, ptol_eff, dtol_eff
         nsteps += 1
     end
 end
@@ -820,3 +897,81 @@ function _unfold_symcover(x::Vector{T}, (u, certified, nsteps, kkt)) where {T}
     xnew = (u[1:n] .- u[n+1:2n]) ./ 2
     return xnew, certified, nsteps, kkt
 end
+
+# Constraints on no cycle of the support graph, found by peeling vertices of
+# degree one. Each is tight with zero multiplier at the minimizer: the vertices
+# beyond it enter no other constraint, so their scales make it tight at no cost
+# and the gradient vanishes there. Seeding them keeps the finish from chasing
+# roundoff-level violations of constraints that are tight by structure.
+function _acyclic_constraints(s::EdgeList, nV::Int)
+    E = length(s.edges)
+    ptr = zeros(Int, nV + 1)
+    for (p, q) in s.edges
+        ptr[p+1] += 1
+        ptr[q+1] += 1
+    end
+    ptr[1] = 1
+    cumsum!(ptr, ptr)
+    adj = zeros(Int, ptr[end] - 1)
+    cursor = ptr[1:nV]
+    for (e, (p, q)) in enumerate(s.edges)
+        adj[cursor[p]] = e
+        cursor[p] += 1
+        adj[cursor[q]] = e
+        cursor[q] += 1
+    end
+    deg = [ptr[v+1] - ptr[v] for v in 1:nV]
+    alive = trues(E)
+    stack = [v for v in 1:nV if deg[v] == 1]
+    out = Int[]
+    while !isempty(stack)
+        v = pop!(stack)
+        deg[v] == 1 || continue
+        for k in ptr[v]:ptr[v+1]-1
+            e = adj[k]
+            alive[e] || continue
+            alive[e] = false
+            push!(out, e)
+            deg[v] -= 1
+            p, q = s.edges[e]
+            w = p == v ? q : p
+            deg[w] -= 1
+            deg[w] == 1 && push!(stack, w)
+        end
+    end
+    return out
+end
+
+function _acyclic_constraints(s::_DiffGridQP, nV::Int)
+    C = s.C
+    L = LinearIndices(C)
+    m, n = size(C)
+    degr = copy(s.nfr)
+    degc = copy(s.nfc)
+    alive = isfinite.(C)
+    out = Int[]
+    rows = [i for i in 1:m if degr[i] == 1]
+    cols = [j for j in 1:n if degc[j] == 1]
+    while !(isempty(rows) && isempty(cols))
+        if !isempty(rows)
+            i = pop!(rows)
+            degr[i] == 1 || continue
+            j = findfirst(view(alive, i, :))
+            alive[i, j] = false
+            push!(out, L[i, j])
+            degr[i] -= 1
+            degc[j] -= 1
+            degc[j] == 1 && push!(cols, j)
+        else
+            j = pop!(cols)
+            degc[j] == 1 || continue
+            i = findfirst(view(alive, :, j))
+            alive[i, j] = false
+            push!(out, L[i, j])
+            degc[j] -= 1
+            degr[i] -= 1
+            degr[i] == 1 && push!(rows, i)
+        end
+    end
+    return out
+end
diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl
index 555e02c..61008ff 100644
--- a/src/minimal_covers.jl
+++ b/src/minimal_covers.jl
@@ -52,10 +52,11 @@ For `Float64`, `:lsqr` uses a Cholesky preconditioner when its predicted storage
 does not exceed `fillbudget` bytes (default `2^30`) and its predicted numeric-
 factorization flop count does not exceed `flopbudget` (default `8e3`) times the
 number of stored entries in `A`'s support, where a symmetric off-diagonal entry
-counts twice, as `nnz` counts it for a fully stored matrix. Otherwise it
-factors a sparser matrix that keeps the couplings only among entries whose
-constraints are active and the diagonal contributions of the rest; whenever
-that factor also exceeds the budgets, it uses a diagonal preconditioner.
+counts twice, as `nnz` counts it for a fully stored matrix. Otherwise, and for
+every other element type, it preconditions with a maximum-weight spanning
+forest of the couplings among entries whose constraints are active, together
+with the diagonal contributions of all entries; that factorization has no fill
+and needs no budget.
 `fillbudget=Inf` and `flopbudget=Inf` together always select the full Cholesky
 preconditioner; `fillbudget=0` or `flopbudget=0` always select the diagonal
 one. The same budgets apply to the active-set method that finishes an
@@ -1382,12 +1383,10 @@ _dense_factor_type(::Type{T}) where {T} = LinearAlgebra.LU{T,Matrix{T},Vector{In
 # the `linsolve` and `precond` choices, the full Cholesky preconditioner's
 # predicted entry count `fill_entries` and flop count `factor_flops` (`0` and
 # `0.0` without a factor), the number of numeric factorizations `nrefactor`
-# (full or heavy-forest), the numbers of LSQR solves
-# preconditioned by the heavy-forest factor (`nforest`) and by the diagonal
-# (`ndiagonal`) when the full factor is over budget, and the number `nzeroed`
-# of positive multipliers zeroed directly on entries with slack while the
-# penalty sat at its cap. `precond` is `:factor` (full factor), `:forest` (the
-# heavy-forest factor served at least one solve), `:diagonal`, or `:none`.
+# (full or forest), the number of LSQR solves preconditioned by the forest
+# factor (`nforest`), and the number `nzeroed` of positive multipliers zeroed
+# directly on entries with slack while the penalty sat at its cap. `precond` is
+# `:factor` (full factor), `:forest`, or `:none`.
 # With `zeroslack=false` the iteration stops at a stall at the penalty cap
 # instead of zeroing such multipliers, for callers that finish the solve by the
 # active-set method of `_polish_cover` or `_polish_symcover`.
@@ -1415,8 +1414,10 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
     use_cg = use_woodbury || use_cg_diff
     ne = supp isa EdgeList ? length(supp.edges) : 0
     use_lsqr = linsolve === :lsqr
-    # CHOLMOD preconditioning is limited to Float64.
-    use_precond = use_lsqr && T === Float64
+    use_precond = use_lsqr
+    # The full-factor preconditioner is a CHOLMOD factorization, so other
+    # element types use the forest preconditioner throughout.
+    use_cholmod = T === Float64
     # Residual multiplicity of each stored entry.
     symmetric = sys.symmetric
     mult = (p, q) -> (symmetric && p != q) ? 2 : 1
@@ -1436,7 +1437,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
     cv = zeros(T, ne + 1)   # √weight · log|A_ij|, with a trailing 0 gauge target
     # Violated entries under the current frozen weights.
     vpat = _violation_pattern(supp)
-    # The dense and factor-preconditioned solves (full or heavy-forest factor)
+    # The dense and factor-preconditioned solves (full or forest factor)
     # refactor only when the weights change, and the weights depend only on
     # (κ, active set): cache that key.
     prevκ = Ref(zero(T))
@@ -1477,7 +1478,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
     # stages. Use its diagonal if the predicted Cholesky factor exceeds the fill
     # or flop budget.
     Msp = _precond_pattern(T, supp, v0, use_precond ? N : 0, mult)
-    MF, fill_entries, flops = use_precond ? _precond_analysis(Msp) : (nothing, 0, 0.0)
+    MF, fill_entries, flops = use_precond && use_cholmod ? _precond_analysis(Msp) : (nothing, 0, 0.0)
     # Number of stored entries of the user's support, counted as `nnz` counts
     # the input matrix: a symmetric off-diagonal entry is stored in both
     # triangles and so counts twice, and every other stored entry counts once.
@@ -1488,7 +1489,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
         end
     end
     flopcap = flopbudget * nstored   # largest affordable predicted flop count
-    use_factor = use_precond && sizeof(T) * fill_entries <= fillbudget &&
+    use_factor = use_precond && use_cholmod && sizeof(T) * fill_entries <= fillbudget &&
                  (flopbudget == Inf || flops <= flopcap)
     # Positions of the entries each factored solve overwrites: the diagonal, and
     # both copies of each off-diagonal support entry.
@@ -1501,18 +1502,18 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
             epos[2 * e] = _nzindex(Msp, q, p)
         end
     end
-    psqrt = zeros(T, use_factor ? 0 : (use_precond ? N : 0))  # `K` of the diagonal preconditioner
-    # Heavy-forest preconditioner, used when the full factor is over budget:
-    # violated entries keep their off-diagonal couplings and the others
-    # contribute only to the diagonal. Its pattern follows the violation
-    # pattern, so each new `(κ, vpat)` gets its own symbolic analysis and budget
-    # check; `forestok` records whether the current one passed.
-    HF = use_precond && !use_factor ? SparseCholesky() : nothing
-    HP = Ref(spzeros(T, 0, 0))
-    forestok = Ref(false)
+    # Forest preconditioner, used when the full factor is over budget or the
+    # element type is not CHOLMOD's: a maximum-weight spanning forest of the
+    # couplings among violated entries, every entry's contribution to the
+    # diagonal, the diagonal of `v0 v0ᵀ`, and a ridge. The couplings among
+    # violated entries form a forest up to a few edges in practice, so little is
+    # lost, and the factorization has no fill, no symbolic analysis, and no
+    # budget. It is rebuilt for each new `(κ, vpat)`.
+    HF = Ref{Union{Nothing,TreeCholesky{T}}}(nothing)
     hI = Int[]
     hJ = Int[]
     hV = T[]
+    hW = T[]
     rhs = zeros(T, use_woodbury ? N : 0, size(U, 2) + 1)
     wsol = similar(rhs)
     WF = use_woodbury ? SparseCholesky() : nothing
@@ -1531,7 +1532,6 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
     nchol = Ref(0)
     nrefactor = Ref(0)
     nforest = Ref(0)
-    ndiagonal = Ref(0)
     solve_weighted = function (x, κ)
         nsolves[] += 1
         if supp isa Grid
@@ -1687,10 +1687,11 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
                 return solve_up!(soly, MF, soly)
             elseif use_precond
                 if prevκ[] != κl || vpat != prevpat
-                    # Heavy-forest matrix: the violated part of `RᵀWR`, the
-                    # diagonal of the rest, the diagonal of `v0 v0ᵀ`, and a
-                    # ridge. `sparse` sums duplicates.
-                    empty!(hI); empty!(hJ); empty!(hV)
+                    # Forest matrix: a maximum-weight spanning forest of the
+                    # violated couplings of `RᵀWR`, the diagonal of all of
+                    # `RᵀWR` and of `v0 v0ᵀ`, and a ridge. The diagonal
+                    # dominates the kept couplings, so the pivots stay positive.
+                    empty!(hI); empty!(hJ); empty!(hV); empty!(hW)
                     fill!(mdiag, zero(T))
                     for (e, (p, q)) in enumerate(edges)
                         w = ws[e]^2
@@ -1700,7 +1701,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
                             mdiag[p] += w
                             mdiag[q] += w
                             if vpat[e]
-                                push!(hI, p, q); push!(hJ, q, p); push!(hV, qs * w, qs * w)
+                                push!(hI, p); push!(hJ, q); push!(hV, qs * w); push!(hW, w)
                             end
                         end
                     end
@@ -1710,27 +1711,15 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
                         dmax = max(dmax, mdiag[p])
                     end
                     ρ = _precond_ridge(dmax)
-                    for p in 1:N
-                        push!(hI, p); push!(hJ, p); push!(hV, mdiag[p] + ρ)
-                    end
-                    HP[] = sparse(hI, hJ, hV, N, N)
-                    HFa = HF::SparseCholesky
-                    analyze!(HFa, HP[])
-                    forestok[] = sizeof(T) * factor_entries(HFa) <= fillbudget &&
-                                 (flopbudget == Inf || factor_flops(HFa) <= flopcap)
-                    if forestok[]
-                        factorize!(HFa, HP[])
-                        nrefactor[] += 1
-                    end
+                    keep = max_weight_forest(N, hI, hJ, hW)
+                    HF[] = tree_cholesky(N, hI[keep], hJ[keep], hV[keep], mdiag .+ ρ)
+                    nrefactor[] += 1
                     prevκ[] = κl
                     copyto!(prevpat, vpat)
                 end
-            end
-            if use_precond && !use_factor && forestok[]
                 nforest[] += 1
-                HFc = HF::SparseCholesky
-                HPc = HP[]
-                # `K = P'L` from the heavy-forest factor, as in the full-factor branch.
+                HFc = HF[]::TreeCholesky{T}
+                # `K = L` from the forest factor, as `K = P'L` in the full-factor branch.
                 Fmul! = function (y, yv)
                     solve_up!(pxv, HFc, yv)
                     for (e, (p, q)) in enumerate(edges)
@@ -1750,55 +1739,11 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
                     solve_ptl!(z, HFc, pg)
                     return z
                 end
-                mul!(px, HPc, x)
-                solve_ptl!(px, HFc, px)
+                lt_mul!(px, HFc, x)
                 soly, it = _lsqr(Fmul!, Ftmul!, cv, px)
                 nlsqr[] += it
                 push!(lsqrtrace, it)
                 return solve_up!(soly, HFc, soly)
-            elseif use_precond && !use_factor
-                # The heavy-forest factor is over budget for this pattern.
-                ndiagonal[] += 1
-                # Weighted degrees, the diagonal of `RᵀWR`.
-                fill!(mdiag, zero(T))
-                for (e, (p, q)) in enumerate(edges)
-                    w = ws[e]^2
-                    if p == q
-                        mdiag[p] += 4 * w
-                    else
-                        mdiag[p] += w
-                        mdiag[q] += w
-                    end
-                end
-                # An unknown outside the support gets an identity row.
-                for p in 1:N
-                    d = mdiag[p] + v0[p]^2
-                    psqrt[p] = sqrt(d > 0 ? d : oneunit(T))
-                end
-                # Diagonal `K` needs only elementwise scaling.
-                Dmul! = function (y, yv)
-                    @. px = yv / psqrt
-                    for (e, (p, q)) in enumerate(edges)
-                        y[e] = ws[e] * (px[p] + qs * px[q])
-                    end
-                    y[g] = dot(v0, px)
-                    return y
-                end
-                Dtmul! = function (z, y)
-                    fill!(pg, zero(T))
-                    for (e, (p, q)) in enumerate(edges)
-                        t = ws[e] * y[e]
-                        pg[p] += t
-                        pg[q] += qs * t
-                    end
-                    @. pg += v0 * y[g]
-                    @. z = pg / psqrt
-                    return z
-                end
-                soly, it = _lsqr(Dmul!, Dtmul!, cv, psqrt .* x)
-                nlsqr[] += it
-                push!(lsqrtrace, it)
-                return soly ./ psqrt
             end
             Amul! = function (y, xx)
                 for (e, (p, q)) in enumerate(edges)
@@ -1966,8 +1911,8 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0;
                drops=Tuple(drops), kkt=Tuple(viols), κs=Tuple(κtrace),
                converged, tol=vtol, polish=nothing,
                linsolve=(use_lsqr ? :lsqr : use_cg ? :woodbury : :dense),
-               precond=(!use_precond ? :none : use_factor ? :factor : nforest[] > 0 ? :forest : :diagonal),
-               nrefactor=nrefactor[], nforest=nforest[], ndiagonal=ndiagonal[], fill_entries, factor_flops=flops, nzeroed)
+               precond=(!use_precond ? :none : use_factor ? :factor : :forest),
+               nrefactor=nrefactor[], nforest=nforest[], fill_entries, factor_flops=flops, nzeroed)
 end
 
 # Worker for `symcover_min(::AbsLog{2})`, returning `(a, stats)`. A supplied
diff --git a/src/stats.jl b/src/stats.jl
index 0f7f420..107f20f 100644
--- a/src/stats.jl
+++ b/src/stats.jl
@@ -102,8 +102,7 @@ otherwise `last(kkt)` and `tol`). `iterations` counts multiplier updates.
 - The remaining fields count the inner linear-algebra work and describe the
   preconditioner; their names and meanings may change between releases:
   `nsolves`, `lsqriters`, `lsqrtrace`, `cgiters`, `cholsolves`, `precond`,
-  `nrefactor`, `nforest`, `ndiagonal`, `fill_entries`, `factor_flops`,
-  `nzeroed`.
+  `nrefactor`, `nforest`, `fill_entries`, `factor_flops`, `nzeroed`.
 """
 struct AugmentedLagrangianStats{T} <: AbstractCoverStats
     converged::Bool
@@ -123,7 +122,6 @@ struct AugmentedLagrangianStats{T} <: AbstractCoverStats
     precond::Symbol
     nrefactor::Int
     nforest::Int
-    ndiagonal::Int
     fill_entries::Int
     factor_flops::Float64
     nzeroed::Int
@@ -183,7 +181,7 @@ function AugmentedLagrangianStats{T}(nt::NamedTuple) where {T}
                                        collect(T, nt.drops), polish, nt.linsolve, nt.nsolves,
                                        nt.lsqriters, collect(Int, nt.lsqrtrace), nt.cgiters,
                                        nt.cholsolves, nt.precond, nt.nrefactor, nt.nforest,
-                                       nt.ndiagonal, nt.fill_entries, nt.factor_flops, nt.nzeroed)
+                                       nt.fill_entries, nt.factor_flops, nt.nzeroed)
 end
 
 LeastSquaresStats(nt::NamedTuple) =
diff --git a/src/tree_cholesky.jl b/src/tree_cholesky.jl
new file mode 100644
index 0000000..234fbbf
--- /dev/null
+++ b/src/tree_cholesky.jl
@@ -0,0 +1,125 @@
+# Cholesky factorization of a symmetric positive definite matrix whose
+# off-diagonal pattern is a forest. Eliminating each tree from its leaves
+# produces no fill, so the factor `L` of `M = L Lᵀ` has one diagonal entry per
+# vertex and one off-diagonal entry per non-root vertex, and factorization,
+# solves, and products are linear in the number of vertices. The arithmetic is
+# that of the element type, so any element type serves.
+struct TreeCholesky{T}
+    parent::Vector{Int}   # parent vertex; 0 for a root
+    order::Vector{Int}    # the vertices, each after all of its children
+    d::Vector{T}          # `L[v, v]`
+    l::Vector{T}          # `L[parent[v], v]`; 0 for a root
+end
+
+# Factor the `N×N` matrix with diagonal `diag` and off-diagonal entries `V[k]`
+# at `(I[k], J[k])` and `(J[k], I[k])`. The edges must form a forest, and the
+# matrix must be positive definite; either failure throws.
+function tree_cholesky(N::Int, I::AbstractVector{Int}, J::AbstractVector{Int},
+                       V::AbstractVector{T}, diag::AbstractVector{T}) where {T}
+    axes(I) == axes(J) == axes(V) ||
+        throw(DimensionMismatch("edge endpoints and values must share axes, got $(axes(I)), $(axes(J)), $(axes(V))"))
+    length(diag) == N || throw(DimensionMismatch("the diagonal has $(length(diag)) entries; expected $N"))
+    ptr = zeros(Int, N + 1)
+    for k in eachindex(I)
+        ptr[I[k]+1] += 1
+        ptr[J[k]+1] += 1
+    end
+    ptr[1] = 1
+    cumsum!(ptr, ptr)
+    adj = zeros(Int, ptr[end] - 1)
+    cursor = ptr[1:N]
+    for k in eachindex(I)
+        adj[cursor[I[k]]] = k
+        cursor[I[k]] += 1
+        adj[cursor[J[k]]] = k
+        cursor[J[k]] += 1
+    end
+    # Breadth-first search from each root lists parents before children.
+    parent = zeros(Int, N)
+    pedge = zeros(Int, N)
+    seen = falses(N)
+    bfs = Int[]
+    sizehint!(bfs, N)
+    for r in 1:N
+        seen[r] && continue
+        seen[r] = true
+        head = length(bfs) + 1
+        push!(bfs, r)
+        while head <= length(bfs)
+            v = bfs[head]
+            head += 1
+            for s in ptr[v]:ptr[v+1]-1
+                k = adj[s]
+                k == pedge[v] && continue
+                w = I[k] == v ? J[k] : I[k]
+                seen[w] && throw(ArgumentError("the edges contain a cycle through vertex $w"))
+                seen[w] = true
+                parent[w] = v
+                pedge[w] = k
+                push!(bfs, w)
+            end
+        end
+    end
+    order = reverse!(bfs)
+    d = zeros(T, N)
+    l = zeros(T, N)
+    s = collect(T, diag)   # each pivot, less the squares of its children's factor entries
+    for v in order
+        s[v] > 0 || throw(ArgumentError("the matrix is not positive definite: pivot $(s[v]) at vertex $v"))
+        d[v] = sqrt(s[v])
+        p = parent[v]
+        if p != 0
+            l[v] = V[pedge[v]] / d[v]
+            s[p] -= l[v]^2
+        end
+    end
+    return TreeCholesky{T}(parent, order, d, l)
+end
+
+# `x = L \ b`; `x` may be `b`.
+function solve_ptl!(x::AbstractVector, F::TreeCholesky, b::AbstractVector)
+    x === b || copyto!(x, b)
+    for v in F.order
+        x[v] /= F.d[v]
+        p = F.parent[v]
+        p == 0 || (x[p] -= F.l[v] * x[v])
+    end
+    return x
+end
+
+# `x = Lᵀ \ b`; `x` may be `b`.
+function solve_up!(x::AbstractVector, F::TreeCholesky, b::AbstractVector)
+    x === b || copyto!(x, b)
+    for v in Iterators.reverse(F.order)
+        p = F.parent[v]
+        p == 0 || (x[v] -= F.l[v] * x[p])
+        x[v] /= F.d[v]
+    end
+    return x
+end
+
+# `y = Lᵀ x`; `y` may be `x`.
+function lt_mul!(y::AbstractVector, F::TreeCholesky, x::AbstractVector)
+    for v in F.order
+        p = F.parent[v]
+        y[v] = F.d[v] * x[v] + (p == 0 ? zero(eltype(y)) : F.l[v] * x[p])
+    end
+    return y
+end
+
+# Indices of the edges of a maximum-weight spanning forest of the graph on
+# `1:N` with edges `(I[k], J[k])` and weights `W[k]`, by Kruskal's algorithm.
+function max_weight_forest(N::Int, I::AbstractVector{Int}, J::AbstractVector{Int}, W::AbstractVector)
+    axes(I) == axes(J) == axes(W) ||
+        throw(DimensionMismatch("edge endpoints and weights must share axes, got $(axes(I)), $(axes(J)), $(axes(W))"))
+    uf = collect(1:N)
+    find(x) = (while uf[x] != x; uf[x] = uf[uf[x]]; x = uf[x]; end; x)
+    keep = Int[]
+    for k in sortperm(W; rev=true)
+        rp, rq = find(I[k]), find(J[k])
+        rp == rq && continue
+        uf[rp] = rq
+        push!(keep, k)
+    end
+    return keep
+end
diff --git a/test/active_set.jl b/test/active_set.jl
index 708de4d..339f673 100644
--- a/test/active_set.jl
+++ b/test/active_set.jl
@@ -159,14 +159,15 @@ end
     b = randn(rng, nV - 1)
     xch = MatrixCovers._laplacian_solve!(zeros(nV - 1), Fch, b)
     xcg = MatrixCovers._laplacian_solve!(zeros(nV - 1), Fcg, b)
-    @test !Fcg.factored[]
+    @test Fcg.fallback[] === nothing
     @test xcg ≈ xch rtol=1e-13
     @test norm(S * xcg - b) <= 8 * eps() * (opnorm(Matrix(S)) * norm(xcg) + norm(b))
+    @test Fcg.err[] <= 8 * eps() * norm(xcg)
     # A solve that CG does not finish falls back to the factorization, which
     # later solves reuse.
     Fcg = MatrixCovers._laplacian_factor(L; flopbudget=0, cgmaxiter=0)
     xfb = MatrixCovers._laplacian_solve!(zeros(nV - 1), Fcg, b)
-    @test Fcg.factored[]
+    @test Fcg.fallback[] !== nothing
     @test xfb ≈ xch rtol=1e-13
     @test MatrixCovers._laplacian_solve!(zeros(nV - 1), Fcg, 2b) ≈ 2xch rtol=1e-13
 end
@@ -190,6 +191,86 @@ end
     @test acg .* bcg' ≈ aref .* bref' rtol=1e-10
 end
 
+@testset "refined conjugate-gradient Laplacian solves" begin
+    # A path Laplacian pinned at one end has condition number of order `nV²`.
+    # CG alone leaves a forward error of that times roundoff; refinement with a
+    # double-double residual brings the solve to the roundoff level of the
+    # solution. The right-hand side is exact, so the error is known.
+    rng = StableRNG(23)
+    function path_laplacian(nV)
+        S = spdiagm(-1 => fill(-1.0, nV - 1), 0 => [fill(2.0, nV - 1); 1.0], 1 => fill(-1.0, nV - 1))
+        return S, sparse(triu(S))
+    end
+    for nV in (300, 3000)
+        S, L = path_laplacian(nV)
+        xtrue = Float64.(rand(rng, -5:5, nV))
+        b = S * xtrue
+        Fcg = MatrixCovers._laplacian_factor(L; flopbudget=0)
+        @test Fcg isa MatrixCovers._LaplacianCG
+        x = MatrixCovers._laplacian_solve!(zeros(nV), Fcg, b)
+        @test Fcg.fallback[] === nothing
+        @test maximum(abs, x - xtrue) <= 64 * eps() * norm(xtrue)
+        @test Fcg.err[] <= 8 * eps() * norm(x)
+        # The residual helper: exact at the solution, accurate for a perturbed one.
+        @test all(iszero, MatrixCovers._residual_dd!(zeros(nV), b, S, xtrue))
+        y = xtrue .+ 1e-9 .* randn(rng, nV)
+        @test MatrixCovers._residual_dd!(zeros(nV), b, S, y) ≈ Vector(S * (xtrue - y)) rtol=1e-5
+    end
+    # Element types CHOLMOD does not factor take the same path, to their own
+    # roundoff, and fall back to a dense Cholesky factorization.
+    nV = 300
+    S, L = path_laplacian(nV)
+    Sb, Lb = BigFloat.(S), BigFloat.(L)
+    xtrue = BigFloat.(rand(rng, -5:5, nV))
+    b = Sb * xtrue
+    Fcg = MatrixCovers._laplacian_factor(Lb)
+    @test Fcg isa MatrixCovers._LaplacianCG{BigFloat}
+    @test Fcg.F === nothing
+    x = MatrixCovers._laplacian_solve!(zeros(BigFloat, nV), Fcg, b)
+    @test Fcg.fallback[] === nothing
+    @test maximum(abs, x - xtrue) <= 64 * eps(BigFloat) * norm(xtrue)
+    @test Fcg.err[] <= 8 * eps(BigFloat) * norm(x)
+    Fcg = MatrixCovers._laplacian_factor(Lb; cgmaxiter=0)
+    x = MatrixCovers._laplacian_solve!(zeros(BigFloat, nV), Fcg, b)
+    @test Fcg.fallback[] isa Cholesky
+    @test Fcg.err[] == 0
+    @test maximum(abs, x - xtrue) <= 1e-6 * norm(xtrue)
+    @test MatrixCovers._laplacian_factor(spzeros(BigFloat, 0, 0)) isa Cholesky
+end
+
+@testset "constraints on no cycle of the support" begin
+    # Rows and columns with a single stored entry, and chains hanging off the
+    # rest of the support, give constraints that are tight with zero multiplier
+    # at the minimizer. The finish activates them before anything else.
+    rng = StableRNG(29)
+    m, n = 14, 12
+    core = [(i, m + j) for j in 1:10 for i in 1:10]
+    # Single-entry rows 11 and 12, a single-entry column 11, and the chain
+    # row 13 – column 12 – row 14 – column 5 into the core.
+    tail = [(11, m + 1), (12, m + 2), (3, m + 11), (13, m + 12), (14, m + 12), (14, m + 5)]
+    edges = vcat(core, tail)
+    cvals = randn(rng, length(edges))
+    nV = m + n
+    s = MatrixCovers.EdgeList{Float64}(edges, cvals, -1.0)
+    acy = MatrixCovers._acyclic_constraints(s, nV)
+    @test sort(acy) == collect(length(core)+1:length(edges))
+    C = fill(-Inf, m, n)
+    for (e, (p, q)) in enumerate(edges)
+        C[p, q-m] = cvals[e]
+    end
+    grid = MatrixCovers._qp_support(MatrixCovers.DiffGrid{Float64}(C, collect(1:m), collect(m+1:m+n)))
+    @test Set(C[MatrixCovers._acyclic_constraints(grid, nV)]) == Set(cvals[length(core)+1:end])
+    # The finish certifies the minimizer of a cover whose support has these
+    # appendages, whether the Laplacians are factored or solved by CG.
+    A = sparse(first.(edges), last.(edges) .- m, exp.(cvals), m, n)
+    aj, bj = MatrixCovers.cover_min_jump(AbsLog{2}(), Matrix(A))
+    for flopbudget in (Inf, 0)
+        a, b, st = _cover_min_abslog2(A; maxouter=2, linsolve=:lsqr, flopbudget)
+        @test st.polish.certified
+        @test a .* b' ≈ aj .* bj' rtol=1e-6
+    end
+end
+
 @testset "difference-grid active-set solver" begin
     # The grid must reproduce the edge list of its supported entries, with rows
     # and columns sharing unknowns and some entries outside the support.
diff --git a/test/element_types.jl b/test/element_types.jl
index 50a256a..6cad9aa 100644
--- a/test/element_types.jl
+++ b/test/element_types.jl
@@ -73,6 +73,20 @@
         @test abig isa Vector{BigFloat}
         @test Float64.(abig) ≈ first(symcover_min(AbsLog{2}(), Float64.(Abig))) rtol=1e-6
         @test abig != BigFloat.(Float32.(abig))
+        # LSQR on such a type is preconditioned by the forest factor, since the
+        # full factor is CHOLMOD's.
+        abig_lsqr, sbig = symcover_min(AbsLog{2}(), sparse(Abig); linsolve=:lsqr)
+        @test MatrixCovers.converged(sbig) && sbig.precond === :forest
+        @test Float64.(abig_lsqr) ≈ Float64.(abig) rtol=1e-6
+        # Double64 goes through the same path at its own precision, and the
+        # forest keeps the iteration count near Float64's.
+        Adbl = Double64.(sparse(Asym))
+        adbl, sdbl = symcover_min(AbsLog{2}(), Adbl; linsolve=:lsqr)
+        a64, s64 = symcover_min(AbsLog{2}(), Float64.(sparse(Asym)); linsolve=:lsqr, fillbudget=0)
+        @test adbl isa Vector{Double64}
+        @test MatrixCovers.converged(sdbl) && sdbl.precond === :forest
+        @test Float64.(adbl) ≈ a64 rtol=1e-8
+        @test sdbl.lsqriters <= 4 * s64.lsqriters
     end
 
     @testset "BigFloat flows through the family" begin
diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl
index 976ebe2..596e7ce 100644
--- a/test/minimal_covers.jl
+++ b/test/minimal_covers.jl
@@ -423,8 +423,8 @@ end
 end
 
 # The LSQR preconditioner runs in two regimes: a Cholesky factor of the weighted
-# normal matrix while the fill budget allows it, and that matrix's diagonal
-# beyond it. Both must reach the same cover.
+# normal matrix while the fill budget allows it, and a spanning forest of the
+# active couplings beyond it. Both must reach the same cover.
 @testset "MMC :lsqr preconditioner regimes" begin
     rng = StableRNG(23)
     n = 200
@@ -437,33 +437,36 @@ end
     af, sf = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr)
     ad, sd = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget=0)
     @test sf.precond === :factor
-    @test sd.precond === :diagonal
-    # The factored regime (re)factorizes the preconditioner at least once and
-    # reports its predicted fill; the diagonal regime never factorizes.
-    @test sf.nrefactor >= 1
+    @test sd.precond === :forest
+    # Both regimes (re)factorize their preconditioner at least once; only the
+    # full factor has predicted fill, and the forest regime still analyzes the
+    # full pattern for the budget check.
+    @test sf.nrefactor >= 1 && sd.nrefactor >= 1
     @test sf.fill_entries > 0
     @test sf.factor_flops > 0
-    @test sd.nrefactor == 0
-    @test sd.factor_flops > 0   # analyzed for the fill check even though the diagonal path is taken
+    @test sd.factor_flops > 0
+    @test sf.nforest == 0 && sd.nforest >= 1
     @test af ≈ ad rtol=1e-6
     @test iscover(af, af, A) && iscover(ad, ad, A)
     @test cover_objective(AbsLog{2}(), af, af, A) ≈ cover_objective(AbsLog{2}(), ad, ad, A) rtol=1e-8
-    # The factored regime is the one that converges in a few iterations per solve.
+    # The full factor is the one that converges in a few iterations per solve.
     @test sf.lsqriters < sd.lsqriters
 
     gf, hf, tf = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr)
     gd, hd, td = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr, fillbudget=0)
-    @test (tf.precond, td.precond) === (:factor, :diagonal)
+    @test (tf.precond, td.precond) === (:factor, :forest)
     @test gf .* hf' ≈ gd .* hd' rtol=1e-6
     @test iscover(gf, hf, G) && iscover(gd, hd, G)
     @test cover_objective(AbsLog{2}(), gf, hf, G) ≈ cover_objective(AbsLog{2}(), gd, hd, G) rtol=1e-8
 
-    # The budget is a keyword of the public solvers, and the paths that never
-    # precondition say so.
+    # The budget is a keyword of the public solvers, the dense path never
+    # preconditions, and an element type CHOLMOD does not factor uses the
+    # forest throughout.
     @test first(symcover_min(AbsLog{2}(), A; linsolve=:lsqr, fillbudget=0)) ≈ ad rtol=1e-6
     @test MatrixCovers._symcover_min_abslog2(A; linsolve=:dense)[2].precond === :none
     Abig = BigFloat.([4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5])
-    @test MatrixCovers._symcover_min_abslog2(Abig; linsolve=:lsqr)[2].precond === :none
+    sbig = MatrixCovers._symcover_min_abslog2(Abig; linsolve=:lsqr)[2]
+    @test sbig.precond === :forest && sbig.fill_entries == 0
 end
 
 """
@@ -501,7 +504,7 @@ end
 # refactorization cost, not just storage: a banded matrix has cheap fill-in and
 # stays on the factor path even at large `n`, while a matrix with the same
 # density but random sparsity fills in enough to cost more to refactorize than
-# the heavy-forest preconditioner it falls back to.
+# the forest preconditioner it falls back to.
 @testset "MMC :lsqr preconditioner cost rule" begin
     rngb = StableRNG(31)
     Aband = banded_sparse_sym(rngb, 2000, 3, 1.0)
@@ -521,19 +524,18 @@ end
     # The in-solver ratio of predicted flops to stored support entries matches
     # `factor_flops / nnz(A)` of the input for a symmetric matrix with both
     # triangles stored: a flop budget just above that ratio selects the factor
-    # path, and just below it selects the heavy-forest path.
+    # path, and just below it selects the forest path.
     ratio = sf.factor_flops / nnz(Arand)
     @test MatrixCovers._symcover_min_abslog2(Arand; linsolve=:lsqr, flopbudget=1.01 * ratio)[2].precond === :factor
     @test MatrixCovers._symcover_min_abslog2(Arand; linsolve=:lsqr, flopbudget=0.99 * ratio)[2].precond === :forest
 
-    @test MatrixCovers._symcover_min_abslog2(Arand; linsolve=:lsqr, flopbudget=0)[2].precond === :diagonal
+    @test MatrixCovers._symcover_min_abslog2(Arand; linsolve=:lsqr, flopbudget=0)[2].precond === :forest
 end
 
-# Random sparse supports have expensive full factors. The heavy-forest
-# preconditioner (violated entries with their couplings, the other entries on
-# the diagonal) must reach the same cover as the full factor, and fall back to
-# the diagonal for patterns whose forest factor is over budget.
-@testset "MMC :lsqr heavy-forest preconditioner" begin
+# Random sparse supports have expensive full factors. The forest
+# preconditioner (a spanning forest of the couplings among violated entries,
+# every entry on the diagonal) must reach the same cover as the full factor.
+@testset "MMC :lsqr forest preconditioner" begin
     rng = MersenneTwister(5)
     n = 400
     S = sprand(rng, n, n, 5 / n) + sparse(randperm(rng, n), 1:n, ones(n), n, n)
@@ -548,7 +550,7 @@ end
         @test tf.precond === :factor
         @test ts.precond === :forest
         @test ts.nforest > 0
-        @test ts.nforest + ts.ndiagonal <= ts.nsolves
+        @test ts.nforest <= ts.nsolves
         @test ts.nrefactor >= 1
         @test gs[Gi] .* hs[Gj] ≈ gf[Gi] .* hf[Gj] rtol=1e-8
         @test iscover(gs, hs, G)
@@ -566,8 +568,8 @@ end
     @test sf.precond === :factor
     @test ss.precond === :forest
     @test ss.nforest > 0
-    @test ss.nforest + ss.ndiagonal <= ss.nsolves
-    @test as[Ai] .* as[Aj] ≈ af[Ai] .* af[Aj] rtol=1e-12
+    @test ss.nforest <= ss.nsolves
+    @test as[Ai] .* as[Aj] ≈ af[Ai] .* af[Aj] rtol=1e-10
     @test iscover(as, as, A)
     @test cover_objective(AbsLog{2}(), as, as, A) ≈ cover_objective(AbsLog{2}(), af, af, A) rtol=1e-12
     cfn = M -> MatrixCovers._symcover_min_abslog2(M; flopbudget=fb)[1]
@@ -595,14 +597,15 @@ end
     end
     @test af[Ai] .* af[Aj] ≈ ga[Ai] .* gb[Aj] rtol=1e-12
 
-    # A flop budget just at what the forest factors need sends some patterns
-    # to the diagonal preconditioner, and the cover is unchanged.
+    # The forest factor needs no budget: any budget below the full factor's
+    # gives the forest regime for every solve, and the cover is unchanged.
     gf, hf, _ = MatrixCovers._cover_min_abslog2(G; flopbudget=Inf)
-    gm, hm, tm = MatrixCovers._cover_min_abslog2(G; flopbudget=1.0)
-    @test tm.nforest > 0 && tm.ndiagonal > 0
-    @test tm.nforest + tm.ndiagonal <= tm.nsolves
-    @test gm[Gi] .* hm[Gj] ≈ gf[Gi] .* hf[Gj] rtol=1e-8
-    @test MatrixCovers._cover_min_abslog2(G; flopbudget=0)[3].nforest == 0
+    for budget in (1.0, 0)
+        gm, hm, tm = MatrixCovers._cover_min_abslog2(G; flopbudget=budget)
+        @test tm.precond === :forest
+        @test 0 < tm.nforest <= tm.nsolves
+        @test gm[Gi] .* hm[Gj] ≈ gf[Gi] .* hf[Gj] rtol=1e-8
+    end
 end
 
 # LSQR continuation starts from the heuristic cover.
@@ -932,7 +935,7 @@ end
     # `2(κ-1)z` per update they hold the KKT residual at 2.6e-6.
     # Without the active-set finish (`polish=false`), slack multipliers are zeroed.
     # Whether the residual then reaches the tolerance depends on rounding (with the
-    # diagonal preconditioner it can stall just above it); the solve warns exactly
+    # forest preconditioner it can stall just above it); the solve warns exactly
     # when it does not. The default solve's finish certifies the result.
     A = banded_sparse_sym(StableRNG(5), 10_000, 3, 1.0)
     for fillbudget in (MatrixCovers.LSQR_FILL_BUDGET, 0)
@@ -941,7 +944,7 @@ end
         end
         @test isempty(logs) == s.converged
         @test all(l -> l.level == Base.CoreLogging.Warn && occursin("try a larger `κ`", l.message), logs)
-        @test s.precond === (fillbudget == 0 ? :diagonal : :factor)
+        @test s.precond === (fillbudget == 0 ? :forest : :factor)
         @test s.nzeroed > 0
         @test maximum(s.κs) == 1e5
         @test s.kkt[end] < 1e-9
diff --git a/test/runtests.jl b/test/runtests.jl
index 68e6a87..96e6e7e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -27,6 +27,7 @@ include("helpers.jl")               # isbalanced, covaries, PENALTIES
     include("powermean.jl")         # the PowerMean soft covers (the soft default)
     include("initializers.jl")      # initialize_symcover/initialize_cover strategies
     include("sparse_cholesky.jl")   # the CHOLMOD driver and the multi-array eachindex helper
+    include("tree_cholesky.jl")     # the fill-free forest factorization behind the LSQR preconditioner
     include("minimal_covers.jl")    # the *_min family (native solvers)
     include("transversal.jl")       # cover_transversal and the maximum-product transversal
     include("active_set.jl")        # exact solution of stalled AbsLog{2} minimal covers
diff --git a/test/tree_cholesky.jl b/test/tree_cholesky.jl
new file mode 100644
index 0000000..f6bbef6
--- /dev/null
+++ b/test/tree_cholesky.jl
@@ -0,0 +1,43 @@
+using MatrixCovers: TreeCholesky, tree_cholesky, max_weight_forest, solve_ptl!, solve_up!, lt_mul!
+
+@testset "tree Cholesky" begin
+    rng = StableRNG(31)
+    for T in (Float64, BigFloat), trial in 1:5
+        N = 40
+        # A random forest: most vertices attach to an earlier one.
+        I = Int[]
+        J = Int[]
+        for v in 2:N
+            rand(rng) < 0.8 && (push!(I, v); push!(J, rand(rng, 1:v-1)))
+        end
+        V = T.(randn(rng, length(I)))
+        # Diagonal dominance makes the matrix positive definite.
+        diag = T.(abs.(randn(rng, N)) .+ 1)
+        for k in eachindex(I)
+            diag[I[k]] += abs(V[k])
+            diag[J[k]] += abs(V[k])
+        end
+        M = Matrix(sparse([I; J; 1:N], [J; I; 1:N], [V; V; diag], N, N))
+        F = tree_cholesky(N, I, J, V, diag)
+        @test F isa TreeCholesky{T}
+        b = T.(randn(rng, N))
+        x = solve_up!(zeros(T, N), F, solve_ptl!(zeros(T, N), F, b))
+        @test M * x ≈ b
+        # `L \ (M x) == Lᵀ x` shows `M == L Lᵀ`.
+        @test solve_ptl!(zeros(T, N), F, M * x) ≈ lt_mul!(zeros(T, N), F, x)
+        # Every operation works in place.
+        xa = copy(b)
+        solve_ptl!(xa, F, xa)
+        solve_up!(xa, F, xa)
+        @test xa ≈ x
+        ya = copy(x)
+        lt_mul!(ya, F, ya)
+        @test ya ≈ lt_mul!(zeros(T, N), F, x)
+    end
+    @test_throws "contain a cycle" tree_cholesky(3, [1, 2, 3], [2, 3, 1], ones(3), fill(10.0, 3))
+    @test_throws "not positive definite" tree_cholesky(2, [1], [2], [5.0], [1.0, 1.0])
+    @test_throws DimensionMismatch tree_cholesky(2, [1], [2], [0.5], [1.0, 1.0, 1.0])
+    # Kruskal keeps the two heaviest edges of the triangle and the pendant edge.
+    @test sort(max_weight_forest(4, [1, 2, 3, 1], [2, 3, 1, 4], [1.0, 3.0, 2.0, 0.5])) == [2, 3, 4]
+    @test max_weight_forest(3, Int[], Int[], Float64[]) == Int[]
+end