diff --git a/Project.toml b/Project.toml
index 104c3cd..e8c67fa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,7 +27,7 @@ Aqua = "0.8"
 ExplicitImports = "1.15"
 HiGHS = "1"
 Ipopt = "1"
-JuMP = "1"
+JuMP = "1.15"
 LinearAlgebra = "1.10.0"
 OffsetArrays = "1"
 PrecompileTools = "1"
diff --git a/README.md b/README.md
index 03f521b..a49ee86 100644
--- a/README.md
+++ b/README.md
@@ -19,5 +19,53 @@ rather than forbid it.  Objective-minimal hard covers (`symcover_min`,
 squared-log-excess penalty is solved natively, with no external solver, while
 the other penalties are available when JuMP and HiGHS (or Ipopt) are loaded.
 
+## Example
+
+```julia
+julia> using MatrixCovers, LinearAlgebra
+
+julia> A = [4.0 2.0; 2.0 16.0];
+
+julia> a = symcover(A)      # a[i] * a[j] >= abs(A[i, j])
+2-element Vector{Float64}:
+ 2.0
+ 4.0
+
+julia> a * a'               # dominates A entrywise
+2×2 Matrix{Float64}:
+ 4.0   8.0
+ 8.0  16.0
+
+julia> iscover(a, A)
+true
+```
+
+Covers are scale-covariant: rescaling the matrix rescales the cover the same way.
+
+```julia
+julia> D = Diagonal([10.0, 0.5]);
+
+julia> symcover(D * A * D) ≈ D * a
+true
+```
+
+Non-symmetric matrices get separate row and column scales from `cover`, and
+`cover_min`/`symcover_min` trade the fast heuristic for a cover that minimizes a
+penalty subject to the same constraint:
+
+```julia
+julia> M = [1.0 2.0 3.0; 6.0 5.0 4.0];
+
+julia> a, b = cover(M);
+
+julia> iscover(a, b, M)
+true
+
+julia> aq, bq = cover_min(AbsLog{2}(), M);   # minimal, solved natively
+
+julia> cover_objective(AbsLog{2}(), aq, bq, M) <= cover_objective(AbsLog{2}(), a, b, M)
+true
+```
+
 See the [documentation](https://HolyLab.github.io/MatrixCovers.jl/dev/)
 for motivation, examples, and a full API reference.
diff --git a/docs/Project.toml b/docs/Project.toml
index a34d103..1bcd2ff 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MatrixCovers = "727e6139-ff52-4636-a344-ed1d23e73ffc"
@@ -9,6 +10,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 Documenter = "1"
 HiGHS = "1"
-JuMP = "1"
+Ipopt = "1"
+JuMP = "1.15"
 LinearAlgebra = "1"
 Unitful = "1"
diff --git a/docs/make.jl b/docs/make.jl
index 5f372da..52b53c2 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,6 +1,6 @@
 using MatrixCovers
 using Documenter
-using JuMP, HiGHS
+using JuMP, HiGHS, Ipopt
 using Unitful   # loaded here so the doctests' own `using` cannot emit precompilation output
 
 DocMeta.setdocmeta!(MatrixCovers, :DocTestSetup, :(using MatrixCovers); recursive=true)
@@ -14,7 +14,10 @@ makedocs(;
         edit_link="main",
         assets=String[],
     ),
-    checkdocs=:exports,
+    # `:public` also covers the bindings marked `public` but not exported (the
+    # extension hooks). It rests on `Base.ispublic`, which Julia 1.10 lacks, so
+    # there the check falls back to the exported set.
+    checkdocs=(VERSION >= v"1.11" ? :public : :exports),
     pages=[
         "Home" => "index.md",
     ],
diff --git a/docs/src/index.md b/docs/src/index.md
index b17da09..6eb8249 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -217,14 +217,30 @@ faster than a general-purpose convex solver.  The other penalties — `AbsLog{1}
 (a linear program) and the non-convex `AbsLinear` variants — are solved through
 [JuMP](https://jump.dev/) and are loaded on demand as a package extension:
 
-```julia
-using JuMP, HiGHS   # HiGHS for AbsLog penalties
-using MatrixCovers
+```jldoctest jumpmin
+julia> using MatrixCovers, JuMP, HiGHS   # HiGHS for the AbsLog penalties
 
-a    = symcover_min(AbsLog{1}(), A)   # L1-minimal symmetric hard cover
-a, b = cover_min(AbsLog{1}(), A)      # L1-minimal general hard cover
+julia> S = [4 1 0; 1 1 5; 0 5 2];      # symmetric
+
+julia> round.(symcover_min(AbsLog{1}(), S); digits=6)   # L1-minimal symmetric hard cover
+3-element Vector{Float64}:
+ 2.0
+ 1.0
+ 5.0
+
+julia> A = [1 2 3; 6 5 4];
+
+julia> a, b = cover_min(AbsLog{1}(), A);   # L1-minimal general hard cover
+
+julia> round.(a * b'; digits=6)            # tight on four of the six entries
+2×3 Matrix{Float64}:
+ 2.4  2.0  3.0
+ 6.0  5.0  7.5
 ```
 
+The solver returns values good to roughly solver tolerance, so these examples
+round before displaying.
+
 The soft `*_min` solvers divide along the same line, but not at the same place:
 [`soft_symcover_min`](@ref) and [`soft_cover_min`](@ref) solve `AbsLog{2}` natively and
 reach for JuMP with Ipopt only for the `AbsLinear` penalties.  They do not accept
@@ -291,15 +307,32 @@ from scratch rather than reading the vector passed in.
 
 For finer control, you can run these manually:
 
-```julia
-using JuMP, Ipopt   # Ipopt for the AbsLinear penalties
-using MatrixCovers
+```jldoctest manualstart
+julia> using MatrixCovers, JuMP, Ipopt   # Ipopt for the AbsLinear penalties
 
-a = symcover_min(AbsLinear{2}(), A)                      # multistart over the whole menu
-a = symcover_min(AbsLinear{2}(), A; strategies=(:geomean,))   # or commit to one start
+julia> S = [4 1 0; 1 1 5; 0 5 2];
 
-a0 = initialize_symcover(A; strategy=:geomean)           # or drive it yourself
-symcover_min!(AbsLinear{2}(), a0, A)
+julia> round.(symcover_min(AbsLinear{2}(), S); digits=6)   # multistart over the whole menu
+3-element Vector{Float64}:
+ 2.0
+ 1.0
+ 5.0
+
+julia> round.(symcover_min(AbsLinear{2}(), S; strategies=(:geomean,)); digits=6)   # or commit to one start
+3-element Vector{Float64}:
+ 2.0
+ 1.0
+ 5.0
+
+julia> a0 = initialize_symcover(S; strategy=:geomean);     # or drive it yourself
+
+julia> symcover_min!(AbsLinear{2}(), a0, S);
+
+julia> round.(a0; digits=6)
+3-element Vector{Float64}:
+ 2.0
+ 1.0
+ 5.0
 ```
 
 The same menu supplies the starting points of the [`soft_symcover`](@ref) and
diff --git a/ext/MatrixCoversIpoptExt.jl b/ext/MatrixCoversIpoptExt.jl
index 5842043..16ebb6f 100644
--- a/ext/MatrixCoversIpoptExt.jl
+++ b/ext/MatrixCoversIpoptExt.jl
@@ -1,15 +1,22 @@
 module MatrixCoversIpoptExt
 
-using JuMP: JuMP, @variable, @objective, @constraint, @NLobjective, @NLconstraint
+using JuMP: JuMP, @variable, @objective, @constraint
 using Ipopt: Ipopt
 using MatrixCovers
 using MatrixCovers: AbsLinear
+using MatrixCovers: _edge_list, _sym_edge_list, _degrees
 
 # The models are built over 1-based positions 1:n; `pr`/`pc` map each position to the
 # corresponding axis index of `A`, and results are scattered back onto vectors
 # whose axes match `A`'s so offset axes are honored. A row/column of `A` with no
 # nonzero entry appears in no constraint or objective term; its scale is set to
 # exactly 0, matching the native solvers.
+#
+# `A` is read through the support hook and gathered into a flat edge list in position
+# space — `ei`/`ej` the endpoints, `elog` the log-magnitude — so a model costs O(nnz)
+# to build rather than O(length(A)). The symmetric list is the full-grid reading; the
+# hard-cover models below sum over its `ei <= ej` half instead, per the objective each
+# one is defined by.
 
 # The AbsLinear objectives are non-convex, so Ipopt returns a local minimum of
 # whichever basin it descends into from the start it is given. That makes the start a
@@ -17,20 +24,35 @@ using MatrixCovers: AbsLinear
 # `*_min!` refiners, which take the start from the caller. The non-mutating `symcover_min`
 # and `cover_min` are multistart drivers over these kernels and live in the main package.
 
-# Ipopt is a local solver on a non-convex problem; any termination other than a
-# solved one means the returned point does not solve the problem posed, so it is an
-# error rather than a result to be quietly handed back.
-function check_solved(model, fname)
-    status = JuMP.termination_status(model)
-    status == JuMP.LOCALLY_SOLVED || status == JuMP.OPTIMAL ||
-        error("$fname: Ipopt terminated with status $status")
-    return nothing
+check_solved(model, fname) =
+    MatrixCovers.check_solved(JuMP.termination_status(model), "Ipopt", fname)
+
+# `set_silent` alone still lets Ipopt print its startup banner, once per session, from
+# its C++ core; `sb` ("suppress banner") is the option that covers it. Every model here
+# is solved for a caller who asked for a cover, not for solver output.
+function _ipopt_model()
+    model = JuMP.Model(Ipopt.Optimizer)
+    JuMP.set_silent(model)
+    JuMP.set_attribute(model, "sb", "yes")
+    return model
+end
+
+# The `i ≤ j` half of a symmetric gather, paired with the multiplicity `w` each entry
+# stands for in the full grid: an off-diagonal pair is two entries of `A`, a diagonal
+# entry one. Carrying the weight is equivalent to summing over both orientations and
+# costs half the terms — and, for the AbsLinear{1} models, half the auxiliary variables.
+function _triangle(fi, fj, flog)
+    keep = [e for e in eachindex(fi) if fi[e] <= fj[e]]
+    return fi[keep], fj[keep], flog[keep], [fi[e] == fj[e] ? 1 : 2 for e in keep]
 end
 
 # ============================================================
 # Hard cover: symcover_min!(::AbsLinear{p}, a, A)
-# Minimizes ∑_{i≤j: A[i,j]≠0} |1 - |A[i,j]|/(a[i]*a[j])|^p
-# subject to a[i]*a[j] ≥ |A[i,j]| for all i≤j with A[i,j]≠0.
+# Minimizes ∑_{i,j: A[i,j]≠0} |1 - |A[i,j]|/(a[i]*a[j])|^p over the full grid, so an
+# off-diagonal pair counts twice and a diagonal entry once — the weighting
+# `cover_objective` reports and the rest of the package minimizes.
+# subject to a[i]*a[j] ≥ |A[i,j]| for all i≤j with A[i,j]≠0, which is the same
+# constraint set as imposing it on the full grid.
 # Variables: α[i] = log(a[i]); constraint: α[i]+α[j] ≥ log|A[i,j]|.
 # ============================================================
 
@@ -40,20 +62,17 @@ function MatrixCovers.symcover_min!(::AbsLinear{2}, a::AbstractVector, A)
     T = float(real(eltype(A)))
     pr = collect(axr)
     n = length(pr)
-    Apos = [A[pr[i], pr[j]] for i in 1:n, j in 1:n]
-    supported = [any(!iszero, @view Apos[i, :]) || any(!iszero, @view Apos[:, i]) for i in 1:n]
-    # Precompute nonzero upper-triangle entries
-    idx = [(i, j) for i in 1:n for j in i:n if Apos[i, j] != 0]
-    logA = Dict((i, j) => log(abs(Apos[i, j])) for (i, j) in idx)
+    fi, fj, flog = _sym_edge_list(A, T)
+    supported = _degrees(fi, n) .> 0
+    ti, tj, tlog, tw = _triangle(fi, fj, flog)
 
-    model = JuMP.Model(Ipopt.Optimizer)
-    JuMP.set_silent(model)
+    model = _ipopt_model()
     start0 = [supported[k] && !iszero(a[pr[k]]) ? log(T(a[pr[k]])) : zero(T) for k in 1:n]
     @variable(model, α[k=1:n], start = start0[k])
-    @NLobjective(model, Min,
-        sum((1 - exp(logA[(i,j)] - α[i] - α[j]))^2 for (i,j) in idx))
-    for (i, j) in idx
-        @constraint(model, α[i] + α[j] >= logA[(i, j)])
+    @objective(model, Min,
+        sum(tw[k] * (1 - exp(tlog[k] - α[ti[k]] - α[tj[k]]))^2 for k in eachindex(ti)))
+    for k in eachindex(ti)
+        @constraint(model, α[ti[k]] + α[tj[k]] >= tlog[k])
     end
     JuMP.optimize!(model)
     check_solved(model, "symcover_min!")
@@ -69,25 +88,22 @@ function MatrixCovers.symcover_min!(::AbsLinear{1}, a::AbstractVector, A)
     T = float(real(eltype(A)))
     pr = collect(axr)
     n = length(pr)
-    Apos = [A[pr[i], pr[j]] for i in 1:n, j in 1:n]
-    supported = [any(!iszero, @view Apos[i, :]) || any(!iszero, @view Apos[:, i]) for i in 1:n]
-    idx = [(i, j) for i in 1:n for j in i:n if Apos[i, j] != 0]
-    logA = Dict((i, j) => log(abs(Apos[i, j])) for (i, j) in idx)
+    fi, fj, flog = _sym_edge_list(A, T)
+    supported = _degrees(fi, n) .> 0
+    ti, tj, tlog, tw = _triangle(fi, fj, flog)
 
-    model = JuMP.Model(Ipopt.Optimizer)
-    JuMP.set_silent(model)
+    model = _ipopt_model()
     start0 = [supported[k] && !iszero(a[pr[k]]) ? log(T(a[pr[k]])) : zero(T) for k in 1:n]
     @variable(model, α[k=1:n], start = start0[k])
     # |1 - exp(lA - αi - αj)| via auxiliary variables t ≥ 0 and slack s
-    @variable(model, t[eachindex(idx)] >= 0)
-    for (k, (i, j)) in enumerate(idx)
-        lA = logA[(i, j)]
-        @NLconstraint(model,  1 - exp(lA - α[i] - α[j]) <= t[k])
-        @NLconstraint(model, -1 + exp(lA - α[i] - α[j]) <= t[k])
+    @variable(model, t[eachindex(ti)] >= 0)
+    for k in eachindex(ti)
+        @constraint(model,  1 - exp(tlog[k] - α[ti[k]] - α[tj[k]]) <= t[k])
+        @constraint(model, -1 + exp(tlog[k] - α[ti[k]] - α[tj[k]]) <= t[k])
     end
-    @objective(model, Min, sum(t))
-    for (i, j) in idx
-        @constraint(model, α[i] + α[j] >= logA[(i, j)])
+    @objective(model, Min, sum(tw[k] * t[k] for k in eachindex(ti)))
+    for k in eachindex(ti)
+        @constraint(model, α[ti[k]] + α[tj[k]] >= tlog[k])
     end
     JuMP.optimize!(model)
     check_solved(model, "symcover_min!")
@@ -114,22 +130,19 @@ function MatrixCovers.cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractV
     T = float(real(eltype(A)))
     pr, pc = collect(axr), collect(axc)
     m, n = length(pr), length(pc)
-    Apos = [A[pr[i], pc[j]] for i in 1:m, j in 1:n]
-    nza = vec(count(!iszero, Apos, dims=2))
-    nzb = vec(count(!iszero, Apos, dims=1))
-    idx = [(i, j) for i in 1:m for j in 1:n if Apos[i, j] != 0]
-    logA = Dict((i, j) => log(abs(Apos[i, j])) for (i, j) in idx)
+    ei, ej, elog = _edge_list(A, T)
+    nza = _degrees(ei, m)
+    nzb = _degrees(ej, n)
 
-    model = JuMP.Model(Ipopt.Optimizer)
-    JuMP.set_silent(model)
+    model = _ipopt_model()
     α0 = [nza[i] > 0 ? log(T(a[pr[i]])) : zero(T) for i in 1:m]
     β0 = [nzb[j] > 0 ? log(T(b[pc[j]])) : zero(T) for j in 1:n]
     @variable(model, α[i=1:m], start = α0[i])
     @variable(model, β[j=1:n], start = β0[j])
-    @NLobjective(model, Min,
-        sum((1 - exp(logA[(i,j)] - α[i] - β[j]))^2 for (i,j) in idx))
-    for (i, j) in idx
-        @constraint(model, α[i] + β[j] >= logA[(i, j)])
+    @objective(model, Min,
+        sum((1 - exp(elog[e] - α[ei[e]] - β[ej[e]]))^2 for e in eachindex(ei)))
+    for e in eachindex(ei)
+        @constraint(model, α[ei[e]] + β[ej[e]] >= elog[e])
     end
     @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n))
     JuMP.optimize!(model)
@@ -149,27 +162,23 @@ function MatrixCovers.cover_min!(::AbsLinear{1}, a::AbstractVector, b::AbstractV
     T = float(real(eltype(A)))
     pr, pc = collect(axr), collect(axc)
     m, n = length(pr), length(pc)
-    Apos = [A[pr[i], pc[j]] for i in 1:m, j in 1:n]
-    nza = vec(count(!iszero, Apos, dims=2))
-    nzb = vec(count(!iszero, Apos, dims=1))
-    idx = [(i, j) for i in 1:m for j in 1:n if Apos[i, j] != 0]
-    logA = Dict((i, j) => log(abs(Apos[i, j])) for (i, j) in idx)
+    ei, ej, elog = _edge_list(A, T)
+    nza = _degrees(ei, m)
+    nzb = _degrees(ej, n)
 
-    model = JuMP.Model(Ipopt.Optimizer)
-    JuMP.set_silent(model)
+    model = _ipopt_model()
     α0 = [nza[i] > 0 ? log(T(a[pr[i]])) : zero(T) for i in 1:m]
     β0 = [nzb[j] > 0 ? log(T(b[pc[j]])) : zero(T) for j in 1:n]
     @variable(model, α[i=1:m], start = α0[i])
     @variable(model, β[j=1:n], start = β0[j])
-    @variable(model, t[eachindex(idx)] >= 0)
-    for (k, (i, j)) in enumerate(idx)
-        lA = logA[(i, j)]
-        @NLconstraint(model,  1 - exp(lA - α[i] - β[j]) <= t[k])
-        @NLconstraint(model, -1 + exp(lA - α[i] - β[j]) <= t[k])
+    @variable(model, t[eachindex(ei)] >= 0)
+    for e in eachindex(ei)
+        @constraint(model,  1 - exp(elog[e] - α[ei[e]] - β[ej[e]]) <= t[e])
+        @constraint(model, -1 + exp(elog[e] - α[ei[e]] - β[ej[e]]) <= t[e])
     end
     @objective(model, Min, sum(t))
-    for (i, j) in idx
-        @constraint(model, α[i] + β[j] >= logA[(i, j)])
+    for e in eachindex(ei)
+        @constraint(model, α[ei[e]] + β[ej[e]] >= elog[e])
     end
     @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n))
     JuMP.optimize!(model)
@@ -196,19 +205,16 @@ function MatrixCovers.soft_symcover_min!(::AbsLinear{2}, a::AbstractVector, A)
     T = float(real(eltype(A)))
     pr = collect(axr)
     n = length(pr)
-    Apos = [A[pr[i], pr[j]] for i in 1:n, j in 1:n]
-    supported = [any(!iszero, @view Apos[i, :]) || any(!iszero, @view Apos[:, i]) for i in 1:n]
-    # Include ALL entries (including zeros, which contribute (1-0)^2=1 regardless of α)
-    nonzero_idx = [(i, j) for i in 1:n, j in 1:n if Apos[i, j] != 0]
-    logA = Dict((i, j) => log(abs(Apos[i, j])) for (i, j) in nonzero_idx)
-    n_zeros = n^2 - length(nonzero_idx)
+    fi, fj, flog = _sym_edge_list(A, T)
+    supported = _degrees(fi, n) .> 0
+    ti, tj, tlog, tw = _triangle(fi, fj, flog)
+    n_zeros = n^2 - length(fi)   # a zero entry contributes (1-0)^2 = 1 regardless of α
 
-    model = JuMP.Model(Ipopt.Optimizer)
-    JuMP.set_silent(model)
+    model = _ipopt_model()
     start0 = [supported[k] ? log(T(a[pr[k]])) : zero(T) for k in 1:n]
     @variable(model, α[k=1:n], start = start0[k])
-    @NLobjective(model, Min,
-        sum((1 - exp(logA[(i,j)] - α[i] - α[j]))^2 for (i,j) in nonzero_idx) + n_zeros)
+    @objective(model, Min,
+        sum(tw[k] * (1 - exp(tlog[k] - α[ti[k]] - α[tj[k]]))^2 for k in eachindex(ti)) + n_zeros)
     JuMP.optimize!(model)
     check_solved(model, "soft_symcover_min!")
     for (i, k) in pairs(pr)
@@ -223,23 +229,20 @@ function MatrixCovers.soft_symcover_min!(::AbsLinear{1}, a::AbstractVector, A)
     T = float(real(eltype(A)))
     pr = collect(axr)
     n = length(pr)
-    Apos = [A[pr[i], pr[j]] for i in 1:n, j in 1:n]
-    supported = [any(!iszero, @view Apos[i, :]) || any(!iszero, @view Apos[:, i]) for i in 1:n]
-    nonzero_idx = [(i, j) for i in 1:n, j in 1:n if Apos[i, j] != 0]
-    logA = Dict((i, j) => log(abs(Apos[i, j])) for (i, j) in nonzero_idx)
-    n_zeros = n^2 - length(nonzero_idx)
+    fi, fj, flog = _sym_edge_list(A, T)
+    supported = _degrees(fi, n) .> 0
+    ti, tj, tlog, tw = _triangle(fi, fj, flog)
+    n_zeros = n^2 - length(fi)   # a zero entry contributes |1 - 0| = 1 regardless of α
 
-    model = JuMP.Model(Ipopt.Optimizer)
-    JuMP.set_silent(model)
+    model = _ipopt_model()
     start0 = [supported[k] ? log(T(a[pr[k]])) : zero(T) for k in 1:n]
     @variable(model, α[k=1:n], start = start0[k])
-    @variable(model, t[eachindex(nonzero_idx)] >= 0)
-    for (k, (i, j)) in enumerate(nonzero_idx)
-        lA = logA[(i, j)]
-        @NLconstraint(model,  1 - exp(lA - α[i] - α[j]) <= t[k])
-        @NLconstraint(model, -1 + exp(lA - α[i] - α[j]) <= t[k])
+    @variable(model, t[eachindex(ti)] >= 0)
+    for k in eachindex(ti)
+        @constraint(model,  1 - exp(tlog[k] - α[ti[k]] - α[tj[k]]) <= t[k])
+        @constraint(model, -1 + exp(tlog[k] - α[ti[k]] - α[tj[k]]) <= t[k])
     end
-    @objective(model, Min, sum(t) + n_zeros)
+    @objective(model, Min, sum(tw[k] * t[k] for k in eachindex(ti)) + n_zeros)
     JuMP.optimize!(model)
     check_solved(model, "soft_symcover_min!")
     for (i, k) in pairs(pr)
@@ -262,21 +265,18 @@ function MatrixCovers.soft_cover_min!(::AbsLinear{2}, a::AbstractVector, b::Abst
     T = float(real(eltype(A)))
     pr, pc = collect(axr), collect(axc)
     m, n = length(pr), length(pc)
-    Apos = [A[pr[i], pc[j]] for i in 1:m, j in 1:n]
-    nza = vec(count(!iszero, Apos, dims=2))
-    nzb = vec(count(!iszero, Apos, dims=1))
-    idx = [(i, j) for i in 1:m for j in 1:n if Apos[i, j] != 0]
-    logA = Dict((i, j) => log(abs(Apos[i, j])) for (i, j) in idx)
-    n_zeros = m * n - length(idx)
+    ei, ej, elog = _edge_list(A, T)
+    nza = _degrees(ei, m)
+    nzb = _degrees(ej, n)
+    n_zeros = m * n - length(ei)
 
-    model = JuMP.Model(Ipopt.Optimizer)
-    JuMP.set_silent(model)
+    model = _ipopt_model()
     α0 = [nza[i] > 0 ? log(T(a[pr[i]])) : zero(T) for i in 1:m]
     β0 = [nzb[j] > 0 ? log(T(b[pc[j]])) : zero(T) for j in 1:n]
     @variable(model, α[i=1:m], start = α0[i])
     @variable(model, β[j=1:n], start = β0[j])
-    @NLobjective(model, Min,
-        sum((1 - exp(logA[(i,j)] - α[i] - β[j]))^2 for (i,j) in idx) + n_zeros)
+    @objective(model, Min,
+        sum((1 - exp(elog[e] - α[ei[e]] - β[ej[e]]))^2 for e in eachindex(ei)) + n_zeros)
     @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n))
     JuMP.optimize!(model)
     check_solved(model, "soft_cover_min!")
@@ -295,24 +295,20 @@ function MatrixCovers.soft_cover_min!(::AbsLinear{1}, a::AbstractVector, b::Abst
     T = float(real(eltype(A)))
     pr, pc = collect(axr), collect(axc)
     m, n = length(pr), length(pc)
-    Apos = [A[pr[i], pc[j]] for i in 1:m, j in 1:n]
-    nza = vec(count(!iszero, Apos, dims=2))
-    nzb = vec(count(!iszero, Apos, dims=1))
-    idx = [(i, j) for i in 1:m for j in 1:n if Apos[i, j] != 0]
-    logA = Dict((i, j) => log(abs(Apos[i, j])) for (i, j) in idx)
-    n_zeros = m * n - length(idx)
+    ei, ej, elog = _edge_list(A, T)
+    nza = _degrees(ei, m)
+    nzb = _degrees(ej, n)
+    n_zeros = m * n - length(ei)
 
-    model = JuMP.Model(Ipopt.Optimizer)
-    JuMP.set_silent(model)
+    model = _ipopt_model()
     α0 = [nza[i] > 0 ? log(T(a[pr[i]])) : zero(T) for i in 1:m]
     β0 = [nzb[j] > 0 ? log(T(b[pc[j]])) : zero(T) for j in 1:n]
     @variable(model, α[i=1:m], start = α0[i])
     @variable(model, β[j=1:n], start = β0[j])
-    @variable(model, t[eachindex(idx)] >= 0)
-    for (k, (i, j)) in enumerate(idx)
-        lA = logA[(i, j)]
-        @NLconstraint(model,  1 - exp(lA - α[i] - β[j]) <= t[k])
-        @NLconstraint(model, -1 + exp(lA - α[i] - β[j]) <= t[k])
+    @variable(model, t[eachindex(ei)] >= 0)
+    for e in eachindex(ei)
+        @constraint(model,  1 - exp(elog[e] - α[ei[e]] - β[ej[e]]) <= t[e])
+        @constraint(model, -1 + exp(elog[e] - α[ei[e]] - β[ej[e]]) <= t[e])
     end
     @objective(model, Min, sum(t) + n_zeros)
     @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n))
diff --git a/ext/MatrixCoversJuMPExt.jl b/ext/MatrixCoversJuMPExt.jl
index b440234..5f1c9b5 100644
--- a/ext/MatrixCoversJuMPExt.jl
+++ b/ext/MatrixCoversJuMPExt.jl
@@ -4,33 +4,43 @@ using JuMP: JuMP, @variable, @objective, @constraint
 using HiGHS: HiGHS
 using MatrixCovers
 using MatrixCovers: AbsLog
+using MatrixCovers: _edge_list, _sym_edge_list, _degrees
 using LinearAlgebra: dot
 
+check_solved(model, fname) =
+    MatrixCovers.check_solved(JuMP.termination_status(model), "HiGHS", fname)
+
 # The models are built over 1-based positions 1:n; `pr`/`pc` map each position to
 # the corresponding axis index of `A`, and results are scattered back onto vectors
 # whose axes match `A`'s so offset axes are honored. A row/column of `A` with no
 # nonzero entry carries no constraint or objective term; its scale is set to exactly
 # 0, matching the native solvers.
+#
+# `A` is read through the support hook and gathered into a flat edge list in position
+# space — `ei`/`ej` the endpoints, `elog` the log-magnitude — so a model costs O(nnz)
+# to build rather than O(length(A)). The symmetric list is the full-grid reading, whose
+# `ei <= ej` half is the constraint set.
 
 # Exact reference for the native `symcover_min(::AbsLog{2})`: same QP, solved by
 # HiGHS. Not exported; used by the test suite to cross-check the native solver.
 function MatrixCovers.symcover_min_jump(::AbsLog{2}, A)
     axr = axes(A, 1)
     axes(A, 2) == axr || throw(ArgumentError("symcover_min_jump requires a square matrix"))
+    MatrixCovers.require_abs_symmetric(A, :symcover_min_jump)
     T = float(real(eltype(A)))
     pr = collect(axr)
     n = length(pr)
-    Apos = [A[pr[i], pr[j]] for i in 1:n, j in 1:n]
-    logA = log.(abs.(Apos))
-    supported = [any(!iszero, @view Apos[i, :]) || any(!iszero, @view Apos[:, i]) for i in 1:n]
+    ei, ej, elog = _sym_edge_list(A, T)
+    supported = _degrees(ei, n) .> 0
     model = JuMP.Model(HiGHS.Optimizer)
     JuMP.set_silent(model)
     @variable(model, α[1:n])
-    @objective(model, Min, sum(abs2, α[i] + α[j] - logA[i, j] for i in 1:n, j in 1:n if Apos[i, j] != 0))
-    for i in 1:n, j in i:n
-        Apos[i, j] != 0 && @constraint(model, α[i] + α[j] - logA[i, j] >= 0)
+    @objective(model, Min, sum(abs2, α[ei[e]] + α[ej[e]] - elog[e] for e in eachindex(ei)))
+    for e in eachindex(ei)
+        ei[e] <= ej[e] && @constraint(model, α[ei[e]] + α[ej[e]] - elog[e] >= 0)
     end
     JuMP.optimize!(model)
+    check_solved(model, "symcover_min_jump")
     a = similar(Array{T}, axr)
     for (i, k) in pairs(pr)
         a[k] = supported[i] ? exp(JuMP.value(α[i])) : zero(T)
@@ -65,13 +75,14 @@ const LEX_L1_SLACK = 1e-9
 # `lin` is the AbsLog{1} objective and `residuals` the expressions α[i]+α[j]-log|A[i,j]| over
 # the support, one per stored entry — the same convention `cover_objective` sums over, so the
 # quadratic minimized here is the AbsLog{2} objective it reports.
-function _minimize_l2_over_l1_face!(model, lin, residuals)
+function _minimize_l2_over_l1_face!(model, lin, residuals, fname)
     isempty(residuals) && return nothing
     linopt = JuMP.value(lin)
     l1 = sum(JuMP.value, residuals)          # the AbsLog{1} objective attained
     @constraint(model, lin <= linopt + LEX_L1_SLACK * max(one(l1), l1))
     @objective(model, Min, sum(r^2 for r in residuals))
     JuMP.optimize!(model)
+    check_solved(model, fname)
     return nothing
 end
 
@@ -83,14 +94,15 @@ end
 function _symcover_min_abslog1(A, start)
     axr = axes(A, 1)
     axes(A, 2) == axr || throw(ArgumentError("symcover_min requires a square matrix"))
+    MatrixCovers.require_abs_symmetric(A, :symcover_min)
     T = float(real(eltype(A)))
     pr = collect(axr)
     n = length(pr)
-    Apos = [A[pr[i], pr[j]] for i in 1:n, j in 1:n]
-    logA = log.(abs.(Apos))
-    colcount = vec(count(!iszero, Apos, dims=1))
-    rowcount = vec(count(!iszero, Apos, dims=2))
-    supported = [colcount[i] > 0 || rowcount[i] > 0 for i in 1:n]
+    ei, ej, elog = _sym_edge_list(A, T)
+    # The gather is already the full grid, so a position's degree is both its row
+    # count and its column count.
+    cnt = _degrees(ei, n)
+    supported = cnt .> 0
     model = JuMP.Model(HiGHS.Optimizer)
     JuMP.set_silent(model)
     if start === nothing
@@ -99,14 +111,15 @@ function _symcover_min_abslog1(A, start)
         α0 = [supported[k] ? log(T(start[pr[k]])) : zero(T) for k in 1:n]
         @variable(model, α[k=1:n], start = α0[k])
     end
-    lin = dot(α, colcount .+ rowcount)
+    lin = dot(α, 2 .* cnt)
     @objective(model, Min, lin)
-    for i in 1:n, j in i:n
-        Apos[i, j] != 0 && @constraint(model, α[i] + α[j] - logA[i, j] >= 0)
+    for e in eachindex(ei)
+        ei[e] <= ej[e] && @constraint(model, α[ei[e]] + α[ej[e]] - elog[e] >= 0)
     end
     JuMP.optimize!(model)
-    residuals = [α[i] + α[j] - logA[i, j] for i in 1:n, j in 1:n if Apos[i, j] != 0]
-    _minimize_l2_over_l1_face!(model, lin, residuals)
+    check_solved(model, "symcover_min")
+    residuals = [α[ei[e]] + α[ej[e]] - elog[e] for e in eachindex(ei)]
+    _minimize_l2_over_l1_face!(model, lin, residuals, "symcover_min")
     a = similar(Array{T}, axr)
     for (i, k) in pairs(pr)
         a[k] = supported[i] ? exp(JuMP.value(α[i])) : zero(T)
@@ -122,19 +135,19 @@ function MatrixCovers.cover_min_jump(::AbsLog{2}, A)
     pc = collect(axc)
     m = length(pr)
     n = length(pc)
-    Apos = [A[pr[i], pc[j]] for i in 1:m, j in 1:n]
-    logA = log.(abs.(Apos))
+    ei, ej, elog = _edge_list(A, T)
     model = JuMP.Model(HiGHS.Optimizer)
     JuMP.set_silent(model)
     @variable(model, α[1:m])
     @variable(model, β[1:n])
-    @objective(model, Min, sum(abs2, α[i] + β[j] - logA[i, j] for i in 1:m, j in 1:n if Apos[i, j] != 0))
-    for i in 1:m, j in 1:n
-        Apos[i, j] != 0 && @constraint(model, α[i] + β[j] - logA[i, j] >= 0)
+    @objective(model, Min, sum(abs2, α[ei[e]] + β[ej[e]] - elog[e] for e in eachindex(ei)))
+    for e in eachindex(ei)
+        @constraint(model, α[ei[e]] + β[ej[e]] - elog[e] >= 0)
     end
-    nza, nzb = vec(sum(!iszero, Apos; dims=2)), vec(sum(!iszero, Apos; dims=1))
+    nza, nzb = _degrees(ei, m), _degrees(ej, n)
     @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n))
     JuMP.optimize!(model)
+    check_solved(model, "cover_min_jump")
     a = similar(Array{T}, axr)
     b = similar(Array{T}, axc)
     for (i, k) in pairs(pr)
@@ -167,10 +180,9 @@ function _cover_min_abslog1(A, start)
     pc = collect(axc)
     m = length(pr)
     n = length(pc)
-    Apos = [A[pr[i], pc[j]] for i in 1:m, j in 1:n]
-    logA = log.(abs.(Apos))
-    colcount = vec(count(!iszero, Apos, dims=1))
-    rowcount = vec(count(!iszero, Apos, dims=2))
+    ei, ej, elog = _edge_list(A, T)
+    rowcount = _degrees(ei, m)
+    colcount = _degrees(ej, n)
     model = JuMP.Model(HiGHS.Optimizer)
     JuMP.set_silent(model)
     if start === nothing
@@ -185,8 +197,8 @@ function _cover_min_abslog1(A, start)
     end
     lin = dot(α, rowcount) + dot(β, colcount)
     @objective(model, Min, lin)
-    for i in 1:m, j in 1:n
-        Apos[i, j] != 0 && @constraint(model, α[i] + β[j] - logA[i, j] >= 0)
+    for e in eachindex(ei)
+        @constraint(model, α[ei[e]] + β[ej[e]] - elog[e] >= 0)
     end
     nza, nzb = rowcount, colcount
     # Gauge pin: the products a[i]*b[j] are unchanged by a -> c*a, b -> b/c, so without this
@@ -194,8 +206,9 @@ function _cover_min_abslog1(A, start)
     # degeneracy the second stage resolves, and stays in force there.
     @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n))
     JuMP.optimize!(model)
-    residuals = [α[i] + β[j] - logA[i, j] for i in 1:m, j in 1:n if Apos[i, j] != 0]
-    _minimize_l2_over_l1_face!(model, lin, residuals)
+    check_solved(model, "cover_min")
+    residuals = [α[ei[e]] + β[ej[e]] - elog[e] for e in eachindex(ei)]
+    _minimize_l2_over_l1_face!(model, lin, residuals, "cover_min")
     a = similar(Array{T}, axr)
     b = similar(Array{T}, axc)
     for (i, k) in pairs(pr)
diff --git a/ext/MatrixCoversSparseArraysExt.jl b/ext/MatrixCoversSparseArraysExt.jl
index e9d9148..3ddaabc 100644
--- a/ext/MatrixCoversSparseArraysExt.jl
+++ b/ext/MatrixCoversSparseArraysExt.jl
@@ -36,9 +36,11 @@ function MatrixCovers.foreach_support_sym(f, A::SparseMatrixCSC)
 end
 
 # Emitted pairs are canonical (row <= col) regardless of uplo: for uplo='L'
-# the stored (i, j) with i >= j is reported as (j, i).
+# the stored (i, j) with i >= j is reported as (j, i). Complex `Hermitian` is
+# admitted alongside the real case because only `abs` of a stored value is ever
+# read, and `abs(A[i,j]) == abs(conj(A[j,i]))`.
 function MatrixCovers.foreach_support_sym(f,
-        S::Union{Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Real,<:SparseMatrixCSC}})
+        S::Union{Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}})
     P = parent(S)
     ax = axes(P, 1)
     axes(P, 2) == ax || throw(DimensionMismatch("foreach_support_sym requires a square matrix, got axes $(axes(P))"))
@@ -61,6 +63,21 @@ function MatrixCovers.foreach_support_sym(f,
     return nothing
 end
 
+# The asymmetric traversal of a wrapped sparse matrix, which the asymmetric cover
+# algorithms and `cover_objective` read even when the matrix is symmetric. Only the
+# named triangle is stored, so each off-diagonal pair is emitted in both
+# orientations and the diagonal once; the magnitudes agree in both, including for a
+# complex `Hermitian`. Without this the wrappers fall back to the generic
+# `AbstractMatrix` method and its full-grid `getindex` scan.
+function MatrixCovers.foreach_support(f,
+        S::Union{Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}})
+    MatrixCovers.foreach_support_sym(S) do i, j, v
+        f(i, j, v)
+        i == j || f(j, i, v)
+    end
+    return nothing
+end
+
 # ============================================================
 # Native minimal-cover (MCM) solvers
 # ============================================================
diff --git a/ext/MatrixCoversUnitfulExt.jl b/ext/MatrixCoversUnitfulExt.jl
index f6d5ef6..88e8877 100644
--- a/ext/MatrixCoversUnitfulExt.jl
+++ b/ext/MatrixCoversUnitfulExt.jl
@@ -10,6 +10,12 @@ const MC = MatrixCovers
 # Exponents of a unit, keyed by atomic unit. Rational, because a symmetric cover
 # halves the diagonal's exponents.
 const UnitExps = Dict{FreeUnits,Rational{Int}}
+
+# A cover objective is dimensionless, so it accumulates in the quantity's own
+# numeric type. This is stated for `Quantity{T}` rather than a concrete quantity
+# type because a matrix whose entries carry different units has exactly that
+# abstract element type.
+MC.scalar_type(::Type{<:Quantity{T}}) where {T} = MC.scalar_type(T)
 const QMatrix = AbstractMatrix{<:Quantity}
 const QVector = AbstractVector{<:Quantity}
 
@@ -17,19 +23,43 @@ const QVector = AbstractVector{<:Quantity}
 # Unit algebra
 # ============================================================
 
+# Unitful exposes no public accessor for the atomic units of a `FreeUnits`, and no
+# supported way to iterate them, so the code below reads its representation
+# directly: `FreeUnits{N,D,A}` carries `N` as a tuple of `Unit{U,D}`, each with
+# fields `tens::Int` and `power::Rational{Int}`. That layout is what the
+# `Unitful.Unit` and `Unitful.FreeUnits` docstrings specify. The decomposition
+# cannot be replaced by unit arithmetic: `gauge` takes a per-atom median over
+# rational exponents, which `*`, `/`, and `^` cannot express.
+
 # An atomic unit at the first power. A prefix belongs to the atom -- `mm` and `m`
 # are distinct -- so a coordinate named in `mm` keeps `mm` in its cover.
 atomic(x::Unit{N,D}) where {N,D} = FreeUnits{(Unit{N,D}(x.tens, 1//1),), D, nothing}()
 
-function exps(u::FreeUnits)
+# The third parameter of `FreeUnits` is the affine offset, `nothing` for an
+# ordinary unit. An affine unit measures from a shifted origin, so `a[i]*b[j]`
+# does not scale it and there is no cover to find; Unitful likewise refuses to
+# multiply affine units.
+function exps(u::FreeUnits{N,D,A}) where {N,D,A}
+    A === nothing || throw(ArgumentError("""
+    affine units are not supported: `$u` measures from a shifted origin, so no \
+    product `a[i]*b[j]` covers an entry in it. Convert to a unit with a true zero \
+    first (`u"K"` for `u"°C"`, `u"Ra"` for `u"°F"`)."""))
     d = UnitExps()
-    for x in typeof(u).parameters[1]
+    for x in N
         a = atomic(x)
         p = get(d, a, 0//1) + x.power
         iszero(p) ? delete!(d, a) : (d[a] = p)
     end
     return d
 end
+
+# `ContextUnits` and `FixedUnits` carry a conversion context this code does not
+# read, so they are refused by name rather than reaching `atomic` as a
+# `MethodError` on an unexported internal.
+exps(u::Unitful.Units) = throw(ArgumentError(
+    "unsupported unit type $(nameof(typeof(u))) for `$u`: MatrixCovers reads `FreeUnits`. " *
+    "Convert with `uconvert(FreeUnits(u), x)`."))
+
 exps(q::Quantity) = exps(unit(q))
 
 # Zero exponents are pruned throughout so that `==` on a `UnitExps` compares
@@ -244,24 +274,24 @@ function asymstart!(f, a::QVector, b::QVector, A::QMatrix, ϕ...; kwargs...)
     return a, b
 end
 
-# Every penalty slot below mirrors MatrixCovers's own method table: where
-# it leaves `ϕ` untyped these do too, and where it dispatches on concrete penalties
-# these enumerate the same ones. Each method is then strictly more specific than the
+# Every penalty slot below mirrors MatrixCovers's own method table: where it
+# accepts any `AbstractCoverPenalty` these do too, and where it dispatches on
+# concrete penalties these enumerate the same ones. Each method is then strictly more specific than the
 # one it shadows -- including those in the JuMP and Ipopt extensions, which leave the
 # matrix slot untyped -- so no ambiguity arises. A penalty the package does not
 # support raises a `MethodError` here exactly as it does on a unitless matrix.
 const PENALTIES = (:(AbsLog{1}), :(AbsLog{2}), :(AbsLinear{1}), :(AbsLinear{2}))
 
-# Heuristic covers and initializers: `ϕ` is untyped upstream and ignored.
+# Heuristic covers and initializers: `ϕ` is checked but not consulted.
 MC.symcover(A::QMatrix; kwargs...) = sym(MC.symcover, A; kwargs...)
-MC.symcover(ϕ, A::QMatrix; kwargs...) = sym(MC.symcover, A, ϕ; kwargs...)
+MC.symcover(ϕ::MC.AbstractCoverPenalty, A::QMatrix; kwargs...) = sym(MC.symcover, A, ϕ; kwargs...)
 MC.symcover!(a::QVector, A::QMatrix; kwargs...) = sym!(MC.symcover!, a, A; kwargs...)
-MC.symcover!(ϕ, a::QVector, A::QMatrix; kwargs...) = sym!(MC.symcover!, a, A, ϕ; kwargs...)
+MC.symcover!(ϕ::MC.AbstractCoverPenalty, a::QVector, A::QMatrix; kwargs...) = sym!(MC.symcover!, a, A, ϕ; kwargs...)
 
 MC.cover(A::QMatrix; kwargs...) = asym(MC.cover, A; kwargs...)
-MC.cover(ϕ, A::QMatrix; kwargs...) = asym(MC.cover, A, ϕ; kwargs...)
+MC.cover(ϕ::MC.AbstractCoverPenalty, A::QMatrix; kwargs...) = asym(MC.cover, A, ϕ; kwargs...)
 MC.cover!(a::QVector, b::QVector, A::QMatrix; kwargs...) = asym!(MC.cover!, a, b, A; kwargs...)
-MC.cover!(ϕ, a::QVector, b::QVector, A::QMatrix; kwargs...) = asym!(MC.cover!, a, b, A, ϕ; kwargs...)
+MC.cover!(ϕ::MC.AbstractCoverPenalty, a::QVector, b::QVector, A::QMatrix; kwargs...) = asym!(MC.cover!, a, b, A, ϕ; kwargs...)
 
 # `cover`/`cover!` dispatch on `Adjoint`/`Transpose` upstream without an eltype
 # bound, so a wrapped `Quantity` matrix needs these to stay unambiguous.
diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl
index de8a318..d4958f7 100644
--- a/src/MatrixCovers.jl
+++ b/src/MatrixCovers.jl
@@ -1,7 +1,7 @@
 module MatrixCovers
 
-using LinearAlgebra: LinearAlgebra, Adjoint, Bidiagonal, Diagonal, SymTridiagonal,
-                     Symmetric, Transpose, Tridiagonal, dot, norm
+using LinearAlgebra: LinearAlgebra, Adjoint, Bidiagonal, Diagonal, Hermitian,
+                     SymTridiagonal, Symmetric, Transpose, Tridiagonal, dot, norm
 using PrecompileTools: PrecompileTools, @compile_workload
 using Random: Random, AbstractRNG, MersenneTwister
 
@@ -15,7 +15,7 @@ export soft_symcover_min, soft_symcover_min!, soft_cover_min, soft_cover_min!
 
 # `public` is parsed as a keyword only from Julia 1.11; this package supports 1.10.
 @static if VERSION >= v"1.11"
-    eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym"))
+    eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type"))
 end
 
 include("penalties.jl")
diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl
index da3fa56..9967fa7 100644
--- a/src/heuristic_covers.jl
+++ b/src/heuristic_covers.jl
@@ -41,7 +41,7 @@ julia> a * a'   # covers |A|: a[i]*a[j] >= abs(A[i, j])
  4.0  4.0
 ```
 """
-symcover(ϕ, A::AbstractMatrix; kwargs...) = symcover(A; kwargs...)
+symcover(ϕ::AbstractCoverPenalty, A::AbstractMatrix; kwargs...) = symcover(A; kwargs...)
 
 function symcover(A::AbstractMatrix; kwargs...)
     axes(A, 2) == axes(A, 1) || throw(ArgumentError("symcover requires a square matrix"))
@@ -61,11 +61,12 @@ must match `axes(A, 1)` (and `A` must be square). `ϕ` has the same meaning as i
 
 See also: [`symcover`](@ref).
 """
-symcover!(ϕ, a::AbstractVector, A::AbstractMatrix; kwargs...) = symcover!(a, A; kwargs...)
+symcover!(ϕ::AbstractCoverPenalty, a::AbstractVector, A::AbstractMatrix; kwargs...) = symcover!(a, A; kwargs...)
 
 function symcover!(a::AbstractVector, A::AbstractMatrix; kwargs...)
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("symcover! requires a square matrix"))
+    require_abs_symmetric(A, :symcover!)
     eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$ax"))
     unconstrained_min!(AbsLog{2}(), a, A)
     boost_feasible!(a, A)
@@ -107,7 +108,7 @@ julia> a * b'
  6.0      5.63709  8.31251
 ```
 """
-cover(ϕ, A::AbstractMatrix; kwargs...) = cover(A; kwargs...)
+cover(ϕ::AbstractCoverPenalty, A::AbstractMatrix; kwargs...) = cover(A; kwargs...)
 
 function cover(A::AbstractMatrix; kwargs...)
     T = float(real(eltype(A)))
@@ -138,7 +139,7 @@ covers.
 
 See also: [`cover`](@ref).
 """
-cover!(ϕ, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) =
+cover!(ϕ::AbstractCoverPenalty, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) =
     cover!(a, b, A; kwargs...)
 
 function cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...)
@@ -237,11 +238,12 @@ function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix
     return nza
 end
 
-function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, b::AbstractVector{T}, A::AbstractMatrix) where T
+function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix)
+    T = float(promote_type(eltype(a), eltype(b)))
     axes(A, 1) == eachindex(a) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, b, A)` requires row indices of `A` to match `a`, got axes(A, 1)=$(axes(A, 1)), axes(a)=$(axes(a))"))
     axes(A, 2) == eachindex(b) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, b, A)` requires column indices of `A` to match `b`, got axes(A, 2)=$(axes(A, 2)), axes(b)=$(axes(b))"))
-    loga = fill!(similar(a), zero(T))
-    logb = fill!(similar(b), zero(T))
+    loga = fill!(similar(a, T), zero(T))
+    logb = fill!(similar(b, T), zero(T))
     nza  = zeros(Int, axes(A, 1))
     nzb  = zeros(Int, axes(A, 2))
     foreach_support(A) do i, j, v
@@ -273,8 +275,13 @@ end
 # every diagonal-zero row is deferred until an off-diagonal neighbor supplies
 # a scale for it (see `boost_feasible_seq!`).
 function init_feasible_diag!(a::AbstractVector{T}, A::AbstractMatrix) where T
-    for k in eachindex(a)
-        a[k] = sqrt(T(abs(A[k, k])))
+    ax = eachindex(a)
+    axes(A) == (ax, ax) || throw(DimensionMismatch("`init_feasible_diag!(a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(axes(A)), axes(a)=$(axes(a)))"))
+    # A diagonal entry the traversal skips is zero, which is the "not yet resolved"
+    # value `boost_feasible_seq!` expects.
+    fill!(a, zero(T))
+    foreach_support_sym(A) do i, j, v
+        i == j && (a[i] = sqrt(T(v)))
     end
     return boost_feasible_seq!(a, A)
 end
@@ -287,10 +294,11 @@ end
 # every entry through it permanently uncoverable, whereas floatmin keeps it
 # representable and, being the smallest normal positive magnitude, changes the
 # resulting cover products negligibly.
-function _tighten_shrink(x::T, lr::T) where T
-    y = x / exp(lr / 2)
+function _tighten_shrink(x, lr)
+    T = float(promote_type(typeof(x), typeof(lr)))
+    y = T(x) / exp(T(lr) / 2)
     if iszero(y) && !iszero(x)
-        y = max(exp(log(x) - lr / 2), floatmin(T))
+        y = max(exp(log(T(x)) - T(lr) / 2), floatmin(T))
     end
     return y
 end
@@ -325,12 +333,13 @@ function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3)
     return a
 end
 
-function tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3) where T
+function tighten_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3)
+    T = float(promote_type(eltype(a), eltype(b)))
     eachindex(a) == axes(A, 1) || throw(DimensionMismatch("indices of a must match row-indexing of A"))
     eachindex(b) == axes(A, 2) || throw(DimensionMismatch("indices of b must match column-indexing of A"))
     lratioa = fill(T(Inf), eachindex(a))
     lratiob = fill(T(Inf), eachindex(b))
-    la, lb = similar(a), similar(b)
+    la, lb = similar(a, T), similar(b, T)
     for _ in 1:maxiter
         map!(log, la, a)   # log(0) = -Inf marks zero scales; see below
         map!(log, lb, b)
@@ -362,11 +371,11 @@ function tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, A::AbstractM
 end
 
 # Adjoint/Transpose wrappers for tighten_cover!.
-function tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, A::Adjoint; kwargs...) where T
+function tighten_cover!(a::AbstractVector, b::AbstractVector, A::Adjoint; kwargs...)
     tighten_cover!(b, a, parent(A); kwargs...)
     return a, b
 end
-function tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, A::Transpose; kwargs...) where T
+function tighten_cover!(a::AbstractVector, b::AbstractVector, A::Transpose; kwargs...)
     tighten_cover!(b, a, parent(A); kwargs...)
     return a, b
 end
@@ -486,7 +495,8 @@ end
 # for every entry visited by `foreach_support`. The
 # diagonal is treated as an ordinary entry. Requires a start with strictly
 # positive scale on every supported row of `a` and column of `b`.
-function boost_feasible!(a::AbstractVector{T}, b::AbstractVector{T}, A::AbstractMatrix) where T
+function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix)
+    T = float(promote_type(eltype(a), eltype(b)))
     IdxA, IdxB = eltype(eachindex(a)), eltype(eachindex(b))
     # `la`/`lb` cache log.(a)/log.(b) and are updated alongside `a`/`b`; see
     # the symmetric method.
@@ -540,30 +550,32 @@ end
 # initialization, a heuristic guaranteed to be O(n^2) seems reasonable.
 function boost_feasible_seq!(a::AbstractVector{T}, A::AbstractMatrix) where T
     ax = eachindex(a)
-    n  = length(ax)
-    f  = first(ax)
-
-    deferred = Tuple{Int,Int,T}[]
-    for j in 1:n-1
-        for ik in 0:n-j-1
-            k   = f + ik
-            l   = k + j
-            Akl = T(abs(A[k, l]))
-            iszero(Akl) && continue
-            ak, al = a[k], a[l]
-            if !iszero(ak) && !iszero(al)
-                aprod = ak * al
-                if aprod < Akl
-                    s = sqrt(Akl / aprod)
-                    a[k] *= s; a[l] *= s
-                end
-            elseif !iszero(ak)
-                a[l] = Akl / ak
-            elseif !iszero(al)
-                a[k] = Akl / al
-            else
-                push!(deferred, (k, l, Akl))
+    axes(A) == (ax, ax) || throw(DimensionMismatch("`boost_feasible_seq!(a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(axes(A)), axes(a)=$(axes(a)))"))
+    I = eltype(ax)
+
+    # The support gathered as off-diagonal pairs, ordered by increasing offset and
+    # then by row, which is the propagation order the heuristic is defined by.
+    pairs = Tuple{I,I,T}[]
+    foreach_support_sym(A) do i, j, v
+        i == j || push!(pairs, (i, j, T(v)))
+    end
+    sort!(pairs; by = ((k, l, _),) -> (l - k, k))
+
+    deferred = Tuple{I,I,T}[]
+    for (k, l, Akl) in pairs
+        ak, al = a[k], a[l]
+        if !iszero(ak) && !iszero(al)
+            aprod = ak * al
+            if aprod < Akl
+                s = sqrt(Akl / aprod)
+                a[k] *= s; a[l] *= s
             end
+        elseif !iszero(ak)
+            a[l] = Akl / ak
+        elseif !iszero(al)
+            a[k] = Akl / al
+        else
+            push!(deferred, (k, l, Akl))
         end
     end
 
@@ -648,7 +660,8 @@ end
 # Requires a start with strictly positive scale on every supported row and column.
 # The shift is covariant under an independent row/column rescaling `D_r*A*D_c`,
 # and is accumulated in the log domain for the reasons given in the symmetric method.
-function inflate_feasible!(a::AbstractVector{T}, b::AbstractVector{T}, A::AbstractMatrix) where T
+function inflate_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix)
+    T = float(promote_type(eltype(a), eltype(b)))
     la, lb = map(log, a), map(log, b)
     tref = Ref(zero(T))
     foreach_support(A) do i, j, v
diff --git a/src/initializers.jl b/src/initializers.jl
index d61ff9d..e56e855 100644
--- a/src/initializers.jl
+++ b/src/initializers.jl
@@ -83,6 +83,7 @@ function initialize_symcover!(a::AbstractVector, A::AbstractMatrix;
                               strategy::Symbol=:hardcover, feasible::Symbol=:inflate, kwargs...)
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("initialize_symcover! requires a square matrix"))
+    require_abs_symmetric(A, :initialize_symcover!)
     eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$ax"))
     _initialize_symcover!(a, A, strategy, feasible; kwargs...) ||
         throw(ArgumentError("strategy=:leaveout requires a support entry that can be dropped without emptying a row"))
@@ -230,20 +231,24 @@ function _leaveout_logmean_init!(a::AbstractVector{T}, A::AbstractMatrix) where
     axes(A) == (ax, ax) || throw(DimensionMismatch("`_leaveout_logmean_init!(a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(axes(A)), axes(a)=$(axes(a)))"))
     nza = unconstrained_min!(AbsLog{2}(), a, A)
     sum(nza) == 0 && return false
+    # One gather serves all three passes below: the two pair scans read the `j >= i`
+    # half of it, the Gauss-Seidel sweeps read whole rows.
+    S = _sym_support(A, T)
     # Most negative residual over the support, with a roundoff-tolerant tie set: exact ties
     # (e.g. z[1,1] == z[2,2] for every 2×2) must not be ordered by floating-point noise.
     zmin = T(Inf)
-    for j in ax, i in first(ax):j
-        Aij = abs(A[i, j])
-        iszero(Aij) && continue
-        zmin = min(zmin, log(T(Aij)) - log(a[i]) - log(a[j]))
+    for i in ax, s in _slots(S, i)
+        j = S.idx[s]
+        j < i && continue
+        zmin = min(zmin, log(S.val[s]) - log(a[i]) - log(a[j]))
     end
     ztol = 64 * eps(T) * max(one(T), abs(zmin))
     ibest = jbest = first(ax) - 1
     Abest = T(Inf)
-    for j in ax, i in first(ax):j
-        Aij = T(abs(A[i, j]))
-        iszero(Aij) && continue
+    for i in ax, s in _slots(S, i)
+        j = S.idx[s]
+        j < i && continue
+        Aij = S.val[s]
         z = log(Aij) - log(a[i]) - log(a[j])
         if z <= zmin + ztol && Aij < Abest
             ibest, jbest, Abest = i, j, Aij
@@ -268,11 +273,10 @@ function _leaveout_logmean_init!(a::AbstractVector{T}, A::AbstractMatrix) where
             iszero(nza[i]) && continue
             num = zero(T)   # Σ_j W[i,j] (log|A[i,j]| - α[j]), α[i]-coefficient split out
             den = zero(T)
-            for j in ax
+            for s in _slots(S, i)
+                j = S.idx[s]
                 (min(i, j) == ibest && max(i, j) == jbest) && continue
-                Aij = abs(A[i, j])
-                iszero(Aij) && continue
-                lAij = log(Aij)
+                lAij = log(S.val[s])
                 if j == i
                     num += lAij
                     den += 2
diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl
index 8f87b16..fcb6401 100644
--- a/src/minimal_covers.jl
+++ b/src/minimal_covers.jl
@@ -203,6 +203,8 @@ end
 function symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_STRATEGIES)
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("symcover_min requires a square matrix"))
+    isempty(strategies) &&
+        throw(ArgumentError("symcover_min: `strategies` must name at least one starting cover"))
     T = float(real(eltype(A)))
     starts = [similar(Array{T}, ax) for _ in strategies]
     # A strategy for which `A` admits no start forfeits its slot; only a menu that yields
@@ -242,9 +244,10 @@ _start_slack(lv::T, li::T, lj::T) where {T} =
 # Shared prologue of the `symcover_min!` kernels: check that the caller's start is a
 # cover of `A`, discard the inert scales on unsupported rows, and move the start onto
 # the coverage boundary exactly, so every kernel begins from a feasible point.
-function _prepare_symcover_start!(a::AbstractVector, A::AbstractMatrix)
+function _prepare_symcover_start!(a::AbstractVector, A::AbstractMatrix, fname=:symcover_min!)
     ax = axes(A, 1)
-    axes(A, 2) == ax || throw(ArgumentError("symcover_min! requires a square matrix"))
+    axes(A, 2) == ax || throw(ArgumentError("$fname requires a square matrix"))
+    require_abs_symmetric(A, fname)
     eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$ax"))
     T = float(eltype(a))
     supp = fill!(similar(a, Bool), false)
@@ -321,7 +324,7 @@ end
 # scalars (‖Mᵀr‖ = ϕbar·α·|c|, ‖r‖ = ϕbar, ‖M‖ from the Frobenius norm of the
 # bidiagonal). Returns `(x, iters)`.
 function _lsqr(Amul!, Atmul!, b::AbstractVector{T}, x0::AbstractVector{T};
-               atol=1e-12, maxiter::Int=2 * (length(b) + length(x0)) + 100) where {T}
+               atol=5000 * eps(T), maxiter::Int=2 * (length(b) + length(x0)) + 100) where {T}
     x = copy(x0)
     u = similar(b)
     Amul!(u, x)
@@ -381,9 +384,13 @@ end
 # unweighted solve; the objective is convex, so it changes the path but not the result.
 function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
                                maxiter::Int=40, linsolve::Symbol=:auto, start=nothing,
-                               boost::Bool=true)
+                               boost::Bool=true, fname=:symcover_min)
     linsolve in (:auto, :dense, :lsqr) ||
         throw(ArgumentError("linsolve must be :auto, :dense, or :lsqr; got :$linsolve"))
+    # The shared entry to the native solve, reached from the sym `*_min` methods in
+    # this package and in the SparseArrays extension, so the precondition is checked
+    # once here rather than at each of them.
+    require_abs_symmetric(A, fname)
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("symcover_min requires a square matrix"))
     # The problem only ever depends on abs.(A), a real quantity, so the working type
@@ -392,27 +399,27 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
     T = float(real(eltype(A)))
     n = length(ax)
     use_lsqr = linsolve === :lsqr
-    # log|A| on the support S; the Newton solve runs on 1-based positions 1:n and
+    # Support entries, one per residual z_ij = α_i + α_j - log|A_ij|, with `cvals`
+    # holding log|A_ij| alongside. The gather reports each off-diagonal pair in both
+    # orientations and the diagonal once, which is the full-grid weighting the
+    # objective is defined with. The Newton solve runs on 1-based positions 1:n and
     # is scattered back onto `a` through `ax` so `A`'s own axes are honored.
-    C = zeros(T, n, n)
-    S = falses(n, n)
-    for (jp, j) in enumerate(ax), (ip, i) in enumerate(ax)
-        Aij = abs(A[i, j])
-        iszero(Aij) && continue
-        C[ip, jp] = log(Aij)
-        S[ip, jp] = true
-    end
-    hassupp = [any(@view S[ip, :]) for ip in 1:n]
-    # Support entries, one per residual z_ij = α_i + α_j - log|A_ij|.
+    G = _sym_support(A, T)
     edges = Tuple{Int,Int}[]
-    for jp in 1:n, ip in 1:n
-        S[ip, jp] && push!(edges, (ip, jp))
+    cvals = T[]
+    hassupp = falses(n)
+    for (ip, i) in enumerate(ax)
+        for s in _slots(G, i)
+            push!(edges, (ip, G.idx[s] - first(ax) + 1))
+            push!(cvals, log(G.val[s]))
+        end
+        hassupp[ip] = !isempty(_slots(G, i))
     end
     ne = length(edges)
     fκ = function (α, κ)
         v = zero(T)
-        for (ip, jp) in edges
-            z = α[ip] + α[jp] - C[ip, jp]
+        for (e, (ip, jp)) in enumerate(edges)
+            z = α[ip] + α[jp] - cvals[e]
             v += (z < 0 ? T(κ) : oneunit(T)) * z^2
         end
         return v
@@ -428,7 +435,6 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
     # of `√W R` (≈ √κ) rather than that of `B` (≈ κ).
     ws = zeros(T, ne)   # √weight per support entry, frozen during one solve
     cv = zeros(T, ne)   # √weight · log|A_ij| (LSQR right-hand side)
-    W = zeros(T, n, n)  # weights (dense path)
     f = zeros(T, n)
     nsolves = Ref(0)
     nlsqr = Ref(0)
@@ -436,10 +442,11 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
         nsolves[] += 1
         if use_lsqr
             for (e, (ip, jp)) in enumerate(edges)
-                w = κ === nothing ? oneunit(T) : ((α[ip] + α[jp] - C[ip, jp]) < 0 ? T(κ) : oneunit(T))
+                c = cvals[e]
+                w = κ === nothing ? oneunit(T) : ((α[ip] + α[jp] - c) < 0 ? T(κ) : oneunit(T))
                 sw = sqrt(w)
                 ws[e] = sw
-                cv[e] = sw * C[ip, jp]
+                cv[e] = sw * c
             end
             Amul! = function (y, x)
                 for (e, (ip, jp)) in enumerate(edges)
@@ -460,16 +467,12 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
             nlsqr[] += it
             return sol
         else
-            fill!(W, zero(T))
             fill!(f, zero(T))
-            for (ip, jp) in edges
-                w = κ === nothing ? oneunit(T) : ((α[ip] + α[jp] - C[ip, jp]) < 0 ? T(κ) : oneunit(T))
-                W[ip, jp] = w
-                f[ip] += w * C[ip, jp]
-            end
             B = zeros(T, n, n)
-            for (ip, jp) in edges
-                w = W[ip, jp]
+            for (e, (ip, jp)) in enumerate(edges)
+                c = cvals[e]
+                w = κ === nothing ? oneunit(T) : ((α[ip] + α[jp] - c) < 0 ? T(κ) : oneunit(T))
+                f[ip] += w * c
                 B[ip, ip] += w
                 B[ip, jp] += w
             end
@@ -494,12 +497,12 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
             αnew = solve_weighted(α, κ)
             t = one(T)
             fnew = fκ(αnew, κ)
-            while fnew > fcur && t > 1e-10
+            while fnew > fcur && t > 500_000 * eps(T)
                 t /= 2
                 fnew = fκ(α .+ t .* (αnew .- α), κ)
             end
             α = α .+ t .* (αnew .- α)
-            fcur - fnew <= 1e-12 * max(fcur, one(T)) && break
+            fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break
             fcur = fnew
         end
     end
@@ -508,9 +511,8 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
     # imposes no coverage constraint and whose optimum the boost would move off.
     γ = zero(T)
     if boost
-        for jp in 1:n, ip in 1:n
-            S[ip, jp] || continue
-            γ = max(γ, (C[ip, jp] - α[ip] - α[jp]) / 2)
+        for (e, (ip, jp)) in enumerate(edges)
+            γ = max(γ, (cvals[e] - α[ip] - α[jp]) / 2)
         end
     end
     # Dense scale vector matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A.
@@ -540,18 +542,27 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
     n = length(axc)
     N = m + n
     use_lsqr = linsolve === :lsqr
-    # log|A| on the support S; internal positions 1:m index rows, m+1:m+n index
-    # columns, and results are scattered back through axr/axc so A's axes are honored.
-    C = zeros(T, m, n)
-    S = falses(m, n)
-    for (jp, j) in enumerate(axc), (ip, i) in enumerate(axr)
-        Aij = abs(A[i, j])
-        iszero(Aij) && continue
-        C[ip, jp] = log(Aij)
-        S[ip, jp] = true
+    # Support entries as edges linking a row position ip to a column position m+jp,
+    # with `cvals` holding log|A_ij| alongside. Internal positions 1:m index rows,
+    # m+1:m+n index columns, and results are scattered back through axr/axc so A's
+    # axes are honored.
+    G = _row_support(A, T)
+    edges = Tuple{Int,Int}[]
+    cvals = T[]
+    nzrow = zeros(Int, m)   # support entries per row, for the balance convention
+    nzcol = zeros(Int, n)   # ditto per column
+    for (ip, i) in enumerate(axr)
+        for s in _slots(G, i)
+            jp = G.idx[s] - first(axc) + 1
+            push!(edges, (ip, m + jp))
+            push!(cvals, log(G.val[s]))
+            nzrow[ip] += 1
+            nzcol[jp] += 1
+        end
     end
-    hasrow = [any(@view S[ip, :]) for ip in 1:m]
-    hascol = [any(@view S[:, jp]) for jp in 1:n]
+    ne = length(edges)
+    hasrow = nzrow .> 0
+    hascol = nzcol .> 0
     # Gauge vector: ±1 on supported variables, 0 on support-free ones (which carry
     # no constraint and are decoupled with an identity row in `solve_weighted`).
     v0 = zeros(T, N)
@@ -561,16 +572,10 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
     for jp in 1:n
         hascol[jp] && (v0[m+jp] = -one(T))
     end
-    # Support entries as edges linking a row position ip to a column position m+jp.
-    edges = Tuple{Int,Int}[]
-    for jp in 1:n, ip in 1:m
-        S[ip, jp] && push!(edges, (ip, m + jp))
-    end
-    ne = length(edges)
     fκ = function (x, κ)
         v = zero(T)
-        for (p, q) in edges
-            z = x[p] + x[q] - C[p, q-m]
+        for (e, (p, q)) in enumerate(edges)
+            z = x[p] + x[q] - cvals[e]
             v += (z < 0 ? T(κ) : oneunit(T)) * z^2
         end
         return v
@@ -583,7 +588,6 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
     # least-squares system so `√W R` has full column rank, applies it matrix-free, and
     # warm-starts from the incoming iterate. After the solve a closed-form shift moves
     # the result to the balance convention, so the pinned gauge is not observable.
-    W = zeros(T, m, n)      # per-support weights (dense path)
     f = zeros(T, N)
     ws = zeros(T, ne)       # √weight per support entry (LSQR path)
     cv = zeros(T, ne + 1)   # √weight · log|A_ij|, with a trailing 0 gauge target
@@ -593,7 +597,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
         nsolves[] += 1
         if use_lsqr
             for (e, (p, q)) in enumerate(edges)
-                c = C[p, q-m]
+                c = cvals[e]
                 w = κ === nothing ? oneunit(T) : ((x[p] + x[q] - c) < 0 ? T(κ) : oneunit(T))
                 sw = sqrt(w)
                 ws[e] = sw
@@ -621,19 +625,13 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
             nlsqr[] += it
             return sol
         else
-            fill!(W, zero(T))
             fill!(f, zero(T))
-            for (p, q) in edges
-                jp = q - m
-                c = C[p, jp]
+            B = v0 * v0'
+            for (e, (p, q)) in enumerate(edges)
+                c = cvals[e]
                 w = κ === nothing ? oneunit(T) : ((x[p] + x[q] - c) < 0 ? T(κ) : oneunit(T))
-                W[p, jp] = w
                 f[p] += w * c
                 f[q] += w * c
-            end
-            B = v0 * v0'
-            for (p, q) in edges
-                w = W[p, q-m]
                 B[p, p] += w
                 B[q, q] += w
                 B[p, q] += w
@@ -681,12 +679,12 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
             xnew = solve_weighted(x, κ)
             t = one(T)
             fnew = fκ(xnew, κ)
-            while fnew > fcur && t > 1e-10
+            while fnew > fcur && t > 500_000 * eps(T)
                 t /= 2
                 fnew = fκ(x .+ t .* (xnew .- x), κ)
             end
             x = x .+ t .* (xnew .- x)
-            fcur - fnew <= 1e-12 * max(fcur, one(T)) && break
+            fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break
             fcur = fnew
         end
     end
@@ -697,25 +695,23 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
     # constraint, and every cover this package returns satisfies it.
     if boost
         γ = zero(T)
-        for jp in 1:n, ip in 1:m
-            S[ip, jp] || continue
-            γ = max(γ, (C[ip, jp] - x[ip] - x[m+jp]) / 2)
+        for (e, (p, q)) in enumerate(edges)
+            γ = max(γ, (cvals[e] - x[p] - x[q]) / 2)
         end
         for p in 1:N
             x[p] += γ
         end
     end
     # Shift along the (e; -e) gauge to the balance convention ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ.
-    nnz = count(S)
     Lα = zero(T)
     Lβ = zero(T)
     for ip in 1:m
-        Lα += count(@view S[ip, :]) * x[ip]
+        Lα += nzrow[ip] * x[ip]
     end
     for jp in 1:n
-        Lβ += count(@view S[:, jp]) * x[m+jp]
+        Lβ += nzcol[jp] * x[m+jp]
     end
-    s = nnz > 0 ? (Lβ - Lα) / (2 * nnz) : zero(T)
+    s = ne > 0 ? (Lβ - Lα) / (2 * ne) : zero(T)
     # Dense scale vectors matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A.
     a = similar(Array{T}, axr)
     b = similar(Array{T}, axc)
@@ -739,7 +735,7 @@ end
 # Both paths inherit the hard workers' handling of a singular signless Laplacian (the
 # `[0 1; 1 0]` support graph among them) and of support-free rows and columns.
 _soft_symcover_min_abslog2(A::AbstractMatrix; kwargs...) =
-    _symcover_min_abslog2(A; κs=(), boost=false, kwargs...)
+    _symcover_min_abslog2(A; κs=(), boost=false, fname=:soft_symcover_min, kwargs...)
 _soft_cover_min_abslog2(A::AbstractMatrix; kwargs...) =
     _cover_min_abslog2(A; κs=(), boost=false, kwargs...)
 
@@ -750,3 +746,17 @@ function symcover_min_jump end
 # Internal exact reference implemented by the MatrixCoversJuMPExt extension; used only to
 # cross-check the native `cover_min(::AbsLog{2})` in the test suite.
 function cover_min_jump end
+
+# A solve that stops for any reason other than a solved one leaves the model holding
+# a point that does not solve the problem posed -- the base of an unbounded ray, or
+# whatever the solver last had. Handing that back would be a minimal cover in name
+# only, so it is an error. `ALMOST_*` statuses are rejected along with the rest:
+# they report a tolerance the caller did not ask for.
+#
+# Takes the status rather than the model so that both solver extensions can share it
+# without the main package depending on JuMP.
+function check_solved(status, solver, fname)
+    Symbol(status) in (:OPTIMAL, :LOCALLY_SOLVED) ||
+        error("$fname: $solver terminated with status $status")
+    return nothing
+end
diff --git a/src/penalties.jl b/src/penalties.jl
index 8286f3a..67d0b47 100644
--- a/src/penalties.jl
+++ b/src/penalties.jl
@@ -73,6 +73,32 @@ struct AbsLinear{p} <: AbstractCoverPenalty end
 # cover_objective
 # ============================================================
 
+"""
+    MatrixCovers.scalar_type(T)
+
+The plain floating-point type underlying the element type `T`, with any units
+removed. [`cover_objective`](@ref) sums the ratios `|A[i,j]| / (a[i]*b[j])`,
+which are dimensionless because a cover requires
+`unit(A[i,j]) == unit(a[i])*unit(b[j])`, so the score is an ordinary number
+whatever the operands carry.
+
+This cannot be expressed as `float(real(T))`: a matrix whose entries carry
+different units has an abstract `eltype`, for which `real` and `oneunit` are
+undefined. A unit-carrying element type therefore needs its own method.
+"""
+scalar_type(::Type{T}) where {T<:Number} = float(real(T))
+
+# The accumulator type for a cover objective over `x`. `eltype` answers it when
+# concrete; a matrix whose entries carry different units can have an element type
+# as abstract as `Quantity`, which names no numeric type at all, so there the
+# elements themselves are consulted. The extra pass is O(length(x)) against an
+# objective that is already O(length(A)).
+function objective_type(x)
+    T = eltype(x)
+    isconcretetype(T) && return scalar_type(T)
+    return mapfoldl(v -> scalar_type(typeof(v)), promote_type, x; init=Bool)
+end
+
 """
     cover_objective(ϕ, a, b, A)
     cover_objective(ϕ, a, A)
@@ -81,29 +107,50 @@ Compute the cover objective `∑_{i,j} ϕ(|A[i,j]| / (a[i] * b[j]))` for the giv
 penalty function `ϕ`. The two-argument form is for symmetric matrices where the cover
 is `a*a'`.
 
+The sum runs over the full grid in both forms, so in the symmetric form each
+off-diagonal pair contributes twice and each diagonal entry once. This weighting
+is what the `sym` solvers minimize, so the score reported here is the quantity
+they optimized; code that reads a symmetric matrix through
+[`foreach_support_sym`](@ref), which reports each pair once, must apply the
+factor of 2 itself to match.
+
 Zero entries of `A` are handled according to `ϕ`:
 - `AbsLog{p}`: zero entries contribute 0 (φ(0) = 0 by convention).
 - `AbsLinear{p}`: zero entries contribute 1 (φ(0) = |1-0|^p = 1).
 
+`eachindex(a)` must match `axes(A, 1)` and `eachindex(b)` must match `axes(A, 2)`.
+
+`A` is read through [`foreach_support`](@ref), so the cost is proportional to the
+support rather than to `length(A)` for a storage type that specializes it.
+
 See also:
 - Penalty types (options for `ϕ`): [`AbsLog`](@ref), [`AbsLinear`](@ref).
 - Solvers: [`symcover`](@ref), [`cover`](@ref), [`soft_symcover`](@ref), [`soft_cover`](@ref).
 """
 function cover_objective(ϕ, a, b, A)
-    T = float(promote_type(eltype(a), eltype(b), real(eltype(A))))
-    s = zero(T)
-    for j in eachindex(b)
-        bj = T(b[j])
-        for i in eachindex(a)
-            ai = T(a[i])
-            Aij = T(abs(A[i, j]))
-            ab = ai * bj
-            # 0/0 → 0 (no cover constraint); nonzero/0 → Inf (violated cover)
-            r = iszero(ab) ? (iszero(Aij) ? zero(T) : typemax(T)) : Aij / ab
-            s += T(ϕ(r))
-        end
+    eachindex(a) == axes(A, 1) ||
+        throw(DimensionMismatch("indices of `a` must match row-indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$(axes(A, 1))"))
+    eachindex(b) == axes(A, 2) ||
+        throw(DimensionMismatch("indices of `b` must match column-indexing of `A`, got eachindex(b)=$(eachindex(b)), axes(A, 2)=$(axes(A, 2))"))
+    T = promote_type(objective_type(A), objective_type(a), objective_type(b))
+    s = Ref(zero(T))
+    nsupport = Ref(0)
+    foreach_support(A) do i, j, v
+        ab = a[i] * b[j]
+        # A nonzero entry over a zero scale is uncovered; `typemax` is the ratio's
+        # stand-in for the infinite excess.
+        r = iszero(ab) ? typemax(T) : T(v / ab)
+        s[] += T(ϕ(r))
+        nsupport[] += 1
     end
-    return s
+    # Every entry outside the support has `r = 0` whatever its scales, including a
+    # zero entry over a zero scale, which constrains nothing. They therefore share
+    # one penalty value, and only their count is needed: zero for `AbsLog`, but a
+    # nonzero constant for `AbsLinear`, which is continuous at `r = 0`.
+    # The guard is not just an optimization: a penalty that is infinite at zero
+    # would otherwise turn a fully-dense `A` into `0 * Inf`, i.e. `NaN`.
+    nzero = length(a) * length(b) - nsupport[]
+    return iszero(nzero) ? s[] : s[] + nzero * T(ϕ(zero(T)))
 end
 cover_objective(ϕ, a, A) = cover_objective(ϕ, a, a, A)
 
diff --git a/src/soft_covers.jl b/src/soft_covers.jl
index 92ab001..f3c72ea 100644
--- a/src/soft_covers.jl
+++ b/src/soft_covers.jl
@@ -74,6 +74,7 @@ soft_symcover(A::AbstractMatrix; kwargs...) = soft_symcover(AbsLinear{2}(), A; k
 soft_symcover(::AbsLog{2}, A::AbstractMatrix; kwargs...) = soft_symcover_min(AbsLog{2}(), A; kwargs...)
 
 function soft_symcover(::AbsLog{1}, A::AbstractMatrix; maxiter::Int=20)
+    require_abs_symmetric(A, :soft_symcover)
     a = soft_symcover_min(AbsLog{2}(), A)   # convex AbsLog{2} minimum: a good start
     _abslog1_iter!(a, A, maxiter)
     return a
@@ -87,12 +88,14 @@ function soft_symcover(::AbsLinear{2}, A::AbstractMatrix; maxiter::Int=32, start
                        rng::AbstractRNG=MersenneTwister(_MULTISTART_SEED))
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("soft_symcover requires a square matrix"))
+    require_abs_symmetric(A, :soft_symcover)
     return _soft_symcover_abslinear2(A, maxiter, starts, _resolve_alias(σ, sigma, 2.0, :σ, :sigma), rng)
 end
 
 function soft_symcover(::AbsLinear{1}, A::AbstractMatrix; maxiter::Int=20, kwargs...)
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("soft_symcover requires a square matrix"))
+    require_abs_symmetric(A, :soft_symcover)
     # Initialize from the AbsLinear{2} soft cover (a good starting point for AbsLinear{1});
     # the multistart is spent there, then the AbsLinear{1} weighted-median descent refines.
     a = soft_symcover(AbsLinear{2}(), A; maxiter=5, kwargs...)
@@ -180,8 +183,8 @@ Supported penalty functions:
 
       u[i] = ∑_j |A[i,j]| v[j] / ∑_j (|A[i,j]| v[j])²   (dually for v[j])
 
-  is monotone, stopping when the relative objective decrease drops below `1e-14` or after
-  `maxiter` sweeps.
+  is monotone, stopping when the relative objective decrease falls to rounding level
+  for the element type, or after `maxiter` sweeps.
 - `AbsLinear{1}()`: initializes from the `AbsLinear{2}()` result, then refines by alternating
   weighted-median updates — each row/column block is minimized exactly, so the descent is
   monotone. Its flat basins are broken by a deterministic lower-median tie-break, giving a
@@ -338,6 +341,8 @@ end
 function soft_symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_STRATEGIES)
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("soft_symcover_min requires a square matrix"))
+    isempty(strategies) &&
+        throw(ArgumentError("soft_symcover_min: `strategies` must name at least one starting cover"))
     T = float(real(eltype(A)))
     starts = [similar(Array{T}, ax) for _ in strategies]
     built = [_initialize_symcover!(a, A, strategy, :none) for (a, strategy) in zip(starts, strategies)]
@@ -387,6 +392,7 @@ end
 function _prepare_soft_symcover_start!(a::AbstractVector, A::AbstractMatrix, fname::Symbol=:soft_symcover_min!)
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("$fname requires a square matrix"))
+    require_abs_symmetric(A, fname)
     eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$ax"))
     supp = fill!(similar(a, Bool), false)
     foreach_support_sym(A) do i, j, v
@@ -554,7 +560,7 @@ end
 # covariance; it must also stay below any genuine basin gap, which is orders of magnitude
 # larger. Candidates stopped short of convergence may differ by far more than roundoff, but
 # such differences are deterministic and co-vary with the frame, so switching on them is safe.
-const _MULTISTART_SWITCHTOL = 1e-9
+_multistart_switchtol(::Type{T}) where {T} = 5_000_000 * eps(T)
 
 # Labeled candidate starts for the symmetric AbsLinear{2} multistart, in selection order.
 # The deterministic starts are the `initialize_symcover` menu: the geometric mean (raw — it is
@@ -585,7 +591,7 @@ function _soft_symcover_abslinear2_inits(A::AbstractMatrix, starts::Int, σ::Rea
         lo = similar(ag)
         _initialize_symcover!(lo, A, :leaveout, :none) && (push!(labels, "leaveout"); push!(inits, lo))
     end
-    if length(inits) < starts && any(iszero, A)
+    if length(inits) < starts && _nsupport(A) < length(A)
         push!(labels, "feasible"); push!(inits, initialize_symcover(A; strategy=:diagfeasible, feasible=:none))
     end
     k = 0
@@ -610,8 +616,9 @@ end
 function _multistart_select(objs)
     besti = firstindex(objs)
     Ebest = objs[besti]
+    switchtol = _multistart_switchtol(float(eltype(objs)))
     for k in eachindex(objs)
-        objs[k] < Ebest * (1 - _MULTISTART_SWITCHTOL) && ((besti, Ebest) = (k, objs[k]))
+        objs[k] < Ebest * (1 - switchtol) && ((besti, Ebest) = (k, objs[k]))
     end
     return besti
 end
@@ -665,18 +672,23 @@ end
 # quantity optimality demands be zero, and it is scale-invariant (each ρ is), so
 # covariant restarts of a rescaled problem exit on the same sweep and the
 # multistart selection stays covariant.
-function _abslinear2_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=1e-8) where T
+function _abslinear2_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=50_000_000 * eps(T)) where T
     ax = eachindex(a)
+    ax == axes(A, 1) || throw(DimensionMismatch("row indices of `A` must match `a`, got $(axes(A, 1)) vs $(ax)"))
+    S = _sym_support(A, T)
     for _ in 1:iter
         maxres = zero(T)
         for k in ax
-            d  = T(abs(A[k, k]))   # diagonal entry
+            d  = zero(T)            # diagonal entry, absent from the support when zero
             s1 = zero(T)            # ∑_{j≠k: A[k,j]≠0} |A[k,j]| / a[j]
             s2 = zero(T)            # ∑_{j≠k: A[k,j]≠0} |A[k,j]|² / a[j]²
-            for j in ax
-                j == k && continue
-                Akj = T(abs(A[k, j]))
-                iszero(Akj) && continue
+            for s in _slots(S, k)
+                j = S.idx[s]
+                Akj = S.val[s]
+                if j == k
+                    d = Akj
+                    continue
+                end
                 inv_aj = one(T) / a[j]
                 s1 += Akj * inv_aj
                 s2 += Akj^2 * inv_aj^2
@@ -749,18 +761,23 @@ end
 # exact fixed point (identical ordering selects the same value), so movement
 # falls to zero there. Relative movement is scale-invariant, so covariant
 # restarts of a rescaled problem exit on the same sweep.
-function _abslinear1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=1e-12) where T
+function _abslinear1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=5000 * eps(T)) where T
     ax  = eachindex(a)
+    ax == axes(A, 1) || throw(DimensionMismatch("row indices of `A` must match `a`, got $(axes(A, 1)) vs $(ax)"))
+    S   = _sym_support(A, T)
     buf = Vector{T}(undef, length(ax))   # reusable buffer for c_j values
     for _ in 1:iter
         maxrel = zero(T)
         for k in ax
-            d  = T(abs(A[k, k]))
+            d  = zero(T)
             nc = 0
-            for j in ax
-                j == k && continue
-                Akj = T(abs(A[k, j]))
-                iszero(Akj) && continue
+            for s in _slots(S, k)
+                j = S.idx[s]
+                Akj = S.val[s]
+                if j == k
+                    d = Akj
+                    continue
+                end
                 aj = a[j]
                 iszero(aj) && continue
                 nc += 1
@@ -803,24 +820,25 @@ end
 # exact fixed point, so movement falls to zero there. Relative movement is
 # scale-invariant, so covariant restarts of a rescaled problem exit on the same
 # sweep.
-function _abslog1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=1e-12) where T
+function _abslog1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=5000 * eps(T)) where T
     ax  = eachindex(a)
+    ax == axes(A, 1) || throw(DimensionMismatch("row indices of `A` must match `a`, got $(axes(A, 1)) vs $(ax)"))
+    S   = _sym_support(A, T)
     buf = Vector{T}(undef, 2 * length(ax) + 1)   # off-diagonals (×1) + diagonal (×2)
     for _ in 1:iter
         maxrel = zero(T)
         for k in ax
             iszero(a[k]) && continue    # zero rows/columns stay uncovered
             n = 0
-            Akk = T(abs(A[k, k]))
-            if !iszero(Akk)
-                lhalf = log(Akk) / 2
-                buf[n += 1] = lhalf
-                buf[n += 1] = lhalf
-            end
-            for j in ax
-                j == k && continue
-                Akj = T(abs(A[k, j]))
-                iszero(Akj) && continue
+            for s in _slots(S, k)
+                j = S.idx[s]
+                Akj = S.val[s]
+                if j == k
+                    lhalf = log(Akj) / 2
+                    buf[n += 1] = lhalf
+                    buf[n += 1] = lhalf
+                    continue
+                end
                 aj = a[j]
                 iszero(aj) && continue
                 buf[n += 1] = log(Akj) - log(aj)
@@ -854,11 +872,14 @@ end
 # `iter` bounds the sweeps; the descent exits early once the largest relative coordinate
 # movement in a sweep drops to `tol` (scale-invariant, so covariant restarts of a rescaled
 # problem exit on the same sweep). Rows/columns with empty support keep scale 0.
-function _abslog1_iter_asym!(a::AbstractVector{T}, b::AbstractVector{T}, A::AbstractMatrix,
-                             iter::Int; tol::Real=1e-12) where T
+function _abslog1_iter_asym!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix,
+                             iter::Int; tol=nothing)
+    T = float(promote_type(eltype(a), eltype(b)))
+    rtol = tol === nothing ? 5000 * eps(T) : T(tol)
     axr, axc = axes(A, 1), axes(A, 2)
     eachindex(a) == axr || throw(DimensionMismatch("row indices of `A` must match `a`, got $(axr) vs $(eachindex(a))"))
     eachindex(b) == axc || throw(DimensionMismatch("column indices of `A` must match `b`, got $(axc) vs $(eachindex(b))"))
+    R, C = _row_support(A, T), _col_support(A, T)
     bufc = Vector{T}(undef, length(axc))   # log-points for an a-row update
     bufr = Vector{T}(undef, length(axr))   # log-points for a b-column update
     for _ in 1:iter
@@ -866,12 +887,10 @@ function _abslog1_iter_asym!(a::AbstractVector{T}, b::AbstractVector{T}, A::Abst
         for i in axr
             iszero(a[i]) && continue       # unsupported rows/columns stay at zero
             nc = 0
-            for j in axc
-                Aij = T(abs(A[i, j]))
-                iszero(Aij) && continue
-                bj = b[j]
+            for s in _slots(R, i)
+                bj = b[R.idx[s]]
                 iszero(bj) && continue
-                bufc[nc += 1] = log(Aij) - log(bj)
+                bufc[nc += 1] = log(R.val[s]) - log(bj)
             end
             nc == 0 && continue
             c = view(bufc, 1:nc)
@@ -885,12 +904,10 @@ function _abslog1_iter_asym!(a::AbstractVector{T}, b::AbstractVector{T}, A::Abst
         for j in axc
             iszero(b[j]) && continue
             nr = 0
-            for i in axr
-                Aij = T(abs(A[i, j]))
-                iszero(Aij) && continue
-                ai = a[i]
+            for s in _slots(C, j)
+                ai = a[C.idx[s]]
                 iszero(ai) && continue
-                bufr[nr += 1] = log(Aij) - log(ai)
+                bufr[nr += 1] = log(C.val[s]) - log(ai)
             end
             nr == 0 && continue
             c = view(bufr, 1:nr)
@@ -901,7 +918,7 @@ function _abslog1_iter_asym!(a::AbstractVector{T}, b::AbstractVector{T}, A::Abst
             iszero(den) || (maxrel = max(maxrel, abs(x - bj) / den))
             b[j] = x
         end
-        maxrel <= T(tol) && break
+        maxrel <= rtol && break
     end
     return a, b
 end
@@ -961,41 +978,42 @@ end
 # exact minimizer. Rows/columns with empty support keep scale 0 and are held fixed.
 # Refines `a`, `b` in place starting from their incoming values.
 #
-# Both half-sweeps accumulate over `i` with `j` held fixed, walking `A` down its columns.
+# Both half-sweeps accumulate over `i` with `j` held fixed, so the support is gathered
+# by column once up front and each sweep walks that gather rather than the full grid.
 #
 # The post-sweep objective costs no extra pass over `A`: once the v-half-sweep has set
 # v[j] = num[j]/den[j] from num[j] = ∑_i M[i,j] u[i] and den[j] = ∑_i (M[i,j] u[i])²,
 # column j contributes
 #     ∑_i (1 - M[i,j] u[i] v[j])² = nnz[j] - 2 v[j] num[j] + v[j]² den[j]
 #                                 = nnz[j] - num[j]²/den[j]
-# with nnz[j] the number of stored nonzeros in column j. An empty column has
-# den[j] = nnz[j] = 0 and contributes nothing.
-#
-# Zero entries of `A` need no branch in the accumulators: M[i,j] = 0 adds zero to both
-# num and den. Only `nnz` distinguishes them, and it is formed once up front.
+# with nnz[j] the number of nonzeros in column j, which is that column's size in the
+# gather. An empty column has den[j] = nnz[j] = 0 and contributes nothing.
 function _mscm_als!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, iter::Int;
-                    tol=1e-14)
+                    tol=nothing)
     axr, axc = axes(A, 1), axes(A, 2)
     eachindex(a) == axr || throw(DimensionMismatch("row indices of `A` must match `a`, got $(axr) vs $(eachindex(a))"))
     eachindex(b) == axc || throw(DimensionMismatch("column indices of `A` must match `b`, got $(axc) vs $(eachindex(b))"))
     T = float(promote_type(eltype(a), eltype(b), real(eltype(A))))
+    # The convergence test is on a relative movement, so its floor is set by the
+    # precision of `T`: a fixed Float64-scaled literal can never be reached in
+    # Float32 (every call would run to `iter`) and stops far short of what a wider
+    # type can resolve.
+    rtol = tol === nothing ? 50 * eps(T) : T(tol)
     # Invert to inverse-scale variables; empty-support rows/columns (scale 0) stay at 0.
     u = map(x -> x > 0 ? inv(T(x)) : zero(T), a)
     v = map(x -> x > 0 ? inv(T(x)) : zero(T), b)
     numu = similar(u)
     denu = similar(u)
-    nnzc = zeros(Int, axc)
-    foreach_support(A) do _, j, _
-        nnzc[j] += 1
-    end
-    E = _mscm_objective(A, u, v)
+    C = _col_support(A, T)
+    E = _mscm_objective(C, u, v)
     for _ in 1:iter
         fill!(numu, zero(T))
         fill!(denu, zero(T))
         for j in axc
             vj = v[j]
-            for i in axr
-                Av = T(abs(A[i, j])) * vj
+            for s in _slots(C, j)
+                i = C.idx[s]
+                Av = C.val[s] * vj
                 numu[i] += Av
                 denu[i] += Av * Av
             end
@@ -1006,17 +1024,17 @@ function _mscm_als!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, ite
         Enew = zero(T)
         for j in axc
             num = den = zero(T)
-            for i in axr
-                Au = T(abs(A[i, j])) * u[i]
+            for s in _slots(C, j)
+                Au = C.val[s] * u[C.idx[s]]
                 num += Au
                 den += Au * Au
             end
             if den > 0
                 v[j] = num / den
-                Enew += nnzc[j] - num * num / den
+                Enew += _ngroup(C, j) - num * num / den
             end
         end
-        E - Enew <= tol * max(E, one(T)) && (E = Enew; break)
+        E - Enew <= rtol * max(E, one(T)) && (E = Enew; break)
         E = Enew
     end
     for i in axr
@@ -1029,15 +1047,14 @@ function _mscm_als!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, ite
 end
 
 # Soft-cover objective in inverse-scale variables: ∑_{i,j: A[i,j]≠0} (1 - |A[i,j]| u[i] v[j])².
-function _mscm_objective(A::AbstractMatrix, u::AbstractVector, v::AbstractVector)
-    T = float(promote_type(real(eltype(A)), eltype(u), eltype(v)))
+# Takes the column-grouped support rather than the matrix, so the sweeps and the
+# objective share one gather.
+function _mscm_objective(C::GroupedSupport{T}, u::AbstractVector, v::AbstractVector) where T
     E = zero(T)
-    for j in axes(A, 2)
+    for j in C.ax
         vj = v[j]
-        for i in axes(A, 1)
-            Aij = T(abs(A[i, j]))
-            iszero(Aij) && continue
-            r = one(T) - Aij * u[i] * vj
+        for s in _slots(C, j)
+            r = one(T) - C.val[s] * u[C.idx[s]] * vj
             E += r * r
         end
     end
@@ -1053,24 +1070,25 @@ end
 # place from their incoming values; exits early once the largest relative coordinate movement
 # in a sweep drops to `tol` (scale-invariant, so covariant restarts of a rescaled problem exit
 # on the same sweep). Rows/columns with empty support keep scale 0.
-function _abslinear1_iter_asym!(a::AbstractVector{T}, b::AbstractVector{T}, A::AbstractMatrix,
-                                iter::Int; tol::Real=1e-12) where T
+function _abslinear1_iter_asym!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix,
+                                iter::Int; tol=nothing)
+    T = float(promote_type(eltype(a), eltype(b)))
+    rtol = tol === nothing ? 5000 * eps(T) : T(tol)
     axr, axc = axes(A, 1), axes(A, 2)
     eachindex(a) == axr || throw(DimensionMismatch("row indices of `A` must match `a`, got $(axr) vs $(eachindex(a))"))
     eachindex(b) == axc || throw(DimensionMismatch("column indices of `A` must match `b`, got $(axc) vs $(eachindex(b))"))
+    R, C = _row_support(A, T), _col_support(A, T)
     bufc = Vector{T}(undef, length(axc))   # c_j buffer for an a-row update
     bufr = Vector{T}(undef, length(axr))   # c_i buffer for a b-column update
     for _ in 1:iter
         maxrel = zero(T)
         for i in axr
             nc = 0
-            for j in axc
-                Aij = T(abs(A[i, j]))
-                iszero(Aij) && continue
-                bj = b[j]
+            for s in _slots(R, i)
+                bj = b[R.idx[s]]
                 iszero(bj) && continue
                 nc += 1
-                bufc[nc] = Aij / bj
+                bufc[nc] = R.val[s] / bj
             end
             x = nc == 0 ? zero(T) : _weighted_self_median!(view(bufc, 1:nc))
             ai  = a[i]
@@ -1080,13 +1098,11 @@ function _abslinear1_iter_asym!(a::AbstractVector{T}, b::AbstractVector{T}, A::A
         end
         for j in axc
             nc = 0
-            for i in axr
-                Aij = T(abs(A[i, j]))
-                iszero(Aij) && continue
-                ai = a[i]
+            for s in _slots(C, j)
+                ai = a[C.idx[s]]
                 iszero(ai) && continue
                 nc += 1
-                bufr[nc] = Aij / ai
+                bufr[nc] = C.val[s] / ai
             end
             x = nc == 0 ? zero(T) : _weighted_self_median!(view(bufr, 1:nc))
             bj  = b[j]
@@ -1094,7 +1110,7 @@ function _abslinear1_iter_asym!(a::AbstractVector{T}, b::AbstractVector{T}, A::A
             iszero(den) || (maxrel = max(maxrel, abs(x - bj) / den))
             b[j] = x
         end
-        maxrel <= T(tol) && break
+        maxrel <= rtol && break
     end
     return a, b
 end
diff --git a/src/support.jl b/src/support.jl
index b436346..0c92486 100644
--- a/src/support.jl
+++ b/src/support.jl
@@ -30,7 +30,8 @@ which must call `f(i, j, abs(A[i, j]))` exactly once for each `(i, j)` with
 `abs(A[i, j]) != 0`, must not call `f` for any other entry (a stored zero is
 still a zero), and must return `nothing`. Emitting an entry twice double-counts
 it in the objective; omitting one silently drops a constraint, yielding a
-"cover" that does not cover.
+"cover" that does not cover. Whatever `f` returns is ignored, so a traversal
+runs to completion and must not be stopped early on the strength of it.
 
 See also: [`foreach_support_sym`](@ref).
 """
@@ -53,18 +54,26 @@ reported in the canonical orientation `i <= j`, the diagonal included, with
 `v = abs(A[i, j])`; zero pairs are skipped. `A` must be square, or a
 `DimensionMismatch` is thrown.
 
-`A` must also be **symmetric in value**, not merely square — this is a
-precondition the function cannot check cheaply and does not try to. It is what
-makes reporting one member of each pair sufficient: a symmetric cover
-constrains `a[i]*a[j]` by a single magnitude, so visiting `(j, i)` as well would
-only duplicate it. Handing this an asymmetric matrix does not error; it silently
-covers a symmetrization of it, and which one is not specified.
+`abs.(A)` must also be **symmetric**, not merely square. That is what makes
+reporting one member of each pair sufficient: a symmetric cover constrains
+`a[i]*a[j]` by a single magnitude, so visiting `(j, i)` as well would only
+duplicate it. Note the predicate is on the magnitudes, so a complex `Hermitian`
+satisfies it — `|A[i,j]| == |conj(A[j,i])|`.
 
-(The `Bidiagonal`/`Tridiagonal` methods, whose two off-diagonal bands are stored
-separately and need not agree, report `max(|A[i,j]|, |A[j,i]|)` — the value a
-full-grid tighten would enforce on both entries. Under the precondition the two
-bands agree and this is just `abs(A[i,j])`; outside it, the choice is
-robustness, not a promise.)
+This traversal does not check the precondition; the public `sym` entry points do,
+before they call it (`MatrixCovers.require_abs_symmetric`).
+
+# Objective weighting
+
+Because each pair is reported once, a caller accumulating a cover objective must
+supply the multiplicity itself: `w = (i == j) ? 1 : 2`. That reproduces the
+`∑_{i,j}` convention of [`cover_objective`](@ref), which runs over the full grid
+and so counts each off-diagonal pair twice and each diagonal entry once. The
+constraint set needs no such correction — `a[i]*a[j] >= |A[i,j]|` and its
+transpose are the same constraint, so imposing it on the `i <= j` triangle alone
+is equivalent to imposing it everywhere. Every solver in this package minimizes
+the full-grid objective, so a cover's reported score and the quantity that was
+minimized agree.
 
 # Extending
 
@@ -74,7 +83,11 @@ To support a new matrix type, define
 
 which must call `f(i, j, v)` exactly once for each pair `i <= j` with
 `v = abs(A[i, j]) != 0`, must not call `f` for zero pairs, and must return
-`nothing`. Reporting the same pair in both orientations double-counts it. A
+`nothing`. Whatever `f` returns is ignored, so a traversal runs to completion
+and must not be stopped early on the strength of it.
+Reporting the same pair in both orientations double-counts it: the off-diagonal
+weight of 2 is the caller's to apply, per *Objective weighting* above, so a pair
+emitted twice is weighted 4. A
 type whose storage is triangular (`Symmetric{<:Any,<:SparseMatrixCSC}` in this
 package's own extension) must map stored `(i, j)` with `i > j` back to `(j, i)`
 rather than emit it as found.
@@ -93,6 +106,165 @@ function foreach_support_sym(f, A::AbstractMatrix)
     return nothing
 end
 
+# Width of the roundoff band the symmetry test allows, in ULPs of the larger of
+# the two magnitudes. A matrix that is symmetric in exact arithmetic can land a
+# ULP or two off after ordinary floating-point work — `D*A*D` for a diagonal
+# rescale is the case that arises throughout this package — and rejecting that
+# would reject input the algorithms handle perfectly well. The band sits far
+# below any genuine asymmetry, which is what the check is for.
+const ASYMMETRY_ULPS = 8
+
+"""
+    MatrixCovers.require_abs_symmetric(A, fname)
+
+Throw unless `abs.(A)` is symmetric to within roundoff, naming `fname` and the
+first offending index pair. Return `nothing` otherwise.
+
+This is the precondition of [`foreach_support_sym`](@ref), enforced at the public
+`sym` entry points rather than inside the traversal, which runs many times per
+solve. An unchecked violation is not a visible failure: the cover returned would
+be a cover of a symmetrization of `A`, plausible-looking and wrong.
+
+The predicate is on the magnitudes rather than on `A` itself, which is both what
+the traversal reads and what admits a complex `Hermitian`.
+"""
+function require_abs_symmetric(A::AbstractMatrix, fname)
+    ax = axes(A, 1)
+    axes(A, 2) == ax ||
+        throw(DimensionMismatch("$fname requires a square matrix, got axes $(axes(A))"))
+    foreach_support(A) do i, j, v
+        w = abs(A[j, i])
+        m = max(v, w)
+        abs(v - w) <= ASYMMETRY_ULPS * eps(float(real(typeof(m)))) * m || throw(ArgumentError("""
+        $fname requires `abs.(A)` to be symmetric, but abs(A[$i,$j]) = $v and \
+        abs(A[$j,$i]) = $w. Wrap `A` in `Symmetric` (or `Hermitian`) to name the \
+        triangle to read; that also skips this check."""))
+    end
+    return nothing
+end
+
+# Storage that makes the precondition structural: the wrapper or the type's own
+# invariant already guarantees `abs(A[i,j]) == abs(A[j,i])`.
+require_abs_symmetric(::Union{Symmetric,Hermitian,Diagonal,SymTridiagonal}, fname) = nothing
+
+# Number of entries `foreach_support` reports, i.e. the size of the stored support.
+function _nsupport(A::AbstractMatrix)
+    n = Ref(0)
+    foreach_support((_, _, _) -> (n[] += 1), A)
+    return n[]
+end
+
+# Support gathered as a flat edge list in 1-based position space, for solver backends
+# whose models are indexed by position rather than by `A`'s own indices: `ei[e]` and
+# `ej[e]` are positions within `axes(A, 1)` and `axes(A, 2)`, and `elog[e]` is
+# `log(abs(A[i,j]))`, the form every model here uses. Building a model from this costs
+# O(nnz) rather than the O(length(A)) of materializing the matrix in position space.
+function _edge_list(A::AbstractMatrix, ::Type{T}) where T
+    or, oc = first(axes(A, 1)) - 1, first(axes(A, 2)) - 1
+    ei, ej, elog = Int[], Int[], T[]
+    foreach_support(A) do i, j, v
+        push!(ei, i - or); push!(ej, j - oc); push!(elog, log(T(v)))
+    end
+    return ei, ej, elog
+end
+
+# Symmetric counterpart: each off-diagonal pair is entered in both orientations and the
+# diagonal once, so the list is the full-grid reading of `A` (see `foreach_support_sym`'s
+# *Objective weighting*). A caller wanting the triangle instead takes the `ei[e] <= ej[e]`
+# half, which is the same constraint set and the other objective convention.
+function _sym_edge_list(A::AbstractMatrix, ::Type{T}) where T
+    o = first(axes(A, 1)) - 1
+    ei, ej, elog = Int[], Int[], T[]
+    foreach_support_sym(A) do i, j, v
+        lv = log(T(v))
+        push!(ei, i - o); push!(ej, j - o); push!(elog, lv)
+        i == j || (push!(ei, j - o); push!(ej, i - o); push!(elog, lv))
+    end
+    return ei, ej, elog
+end
+
+# Number of edges incident on each of `n` positions, counting one endpoint per edge.
+function _degrees(endpoints, n)
+    d = zeros(Int, n)
+    for p in endpoints
+        d[p] += 1
+    end
+    return d
+end
+
+# Support of `A` gathered into per-group neighbor lists, in compressed form: the
+# entries of group `g` occupy the slots `_slots(S, g)`, with `S.idx[s]` the
+# partner index and `S.val[s]` the magnitude.
+#
+# A coordinate-descent kernel revisits every row on every sweep. Reading them
+# through `foreach_support` directly would mean one traversal per sweep, and a
+# traversal cannot be restarted mid-row; gathering once up front costs O(nnz)
+# storage and turns each sweep into a walk over the support instead of over the
+# full grid.
+#
+# `ptr` is indexed by position within `ax` rather than by the index itself, so
+# offset axes need no special case.
+struct GroupedSupport{T,I,R<:AbstractUnitRange}
+    ax::R
+    ptr::Vector{Int}
+    idx::Vector{I}
+    val::Vector{T}
+end
+
+_slots(S::GroupedSupport, g) = (p = g - first(S.ax) + 1; S.ptr[p]:S.ptr[p+1]-1)
+
+# Number of stored entries in group `g`.
+_ngroup(S::GroupedSupport, g) = length(_slots(S, g))
+
+# Build from a traversal. `each` is called twice with a callback `g(group,
+# partner, v)`: once to count each group's entries, once to fill them.
+function _grouped_support(each, ax::AbstractUnitRange, ::Type{I}, ::Type{T}) where {I,T}
+    off = first(ax) - 1
+    ptr = zeros(Int, length(ax) + 1)
+    each() do g, _, _
+        ptr[g-off+1] += 1
+    end
+    ptr[1] = 1
+    cumsum!(ptr, ptr)
+    n = ptr[end] - 1
+    idx = Vector{I}(undef, n)
+    val = Vector{T}(undef, n)
+    cursor = ptr[1:end-1]        # next free slot of each group, by position
+    each() do g, p, v
+        s = cursor[g-off]
+        idx[s] = p
+        val[s] = T(v)
+        cursor[g-off] = s + 1
+    end
+    return GroupedSupport(ax, ptr, idx, val)
+end
+
+# Group the support of `A` by row; partners are column indices.
+_row_support(A::AbstractMatrix, ::Type{T}) where {T} =
+    _grouped_support(axes(A, 1), eltype(axes(A, 2)), T) do g
+        foreach_support((i, j, v) -> g(i, j, v), A)
+    end
+
+# Group the support of `A` by column; partners are row indices.
+_col_support(A::AbstractMatrix, ::Type{T}) where {T} =
+    _grouped_support(axes(A, 2), eltype(axes(A, 1)), T) do g
+        foreach_support((i, j, v) -> g(j, i, v), A)
+    end
+
+# Group the symmetric support of `A` by row: group `i` holds every `j` with
+# `abs(A[i,j]) != 0`, the off-diagonal pairs entered in both orientations and the
+# diagonal once. That is the full-grid reading of a row, so a kernel accumulating
+# over these groups gets the `∑_{i,j}` weighting of `cover_objective` — each
+# off-diagonal pair twice, each diagonal entry once — without applying a
+# multiplicity factor of its own.
+_sym_support(A::AbstractMatrix, ::Type{T}) where {T} =
+    _grouped_support(axes(A, 1), eltype(axes(A, 1)), T) do g
+        foreach_support_sym(A) do i, j, v
+            g(i, j, v)
+            i == j || g(j, i, v)
+        end
+    end
+
 function foreach_support(f, D::Diagonal)
     for i in axes(D, 1)
         v = abs(D[i, i])
@@ -163,9 +335,8 @@ function foreach_support(f, A::Tridiagonal)
     return nothing
 end
 
-# Bidiagonal/Tridiagonal symmetric-contract traversal: a structural-zero band
-# entry (e.g. the subdiagonal of an upper Bidiagonal) reads as exact zero via
-# getindex, so `max` over both entries is correct without special-casing uplo.
+# The two band entries agree to within `ASYMMETRY_ULPS`, so `max` reports the one
+# that dominates both without special-casing uplo.
 function foreach_support_sym(f, A::Union{Bidiagonal,Tridiagonal})
     ax = axes(A, 1)
     for i in ax
diff --git a/test/element_types.jl b/test/element_types.jl
new file mode 100644
index 0000000..790be06
--- /dev/null
+++ b/test/element_types.jl
@@ -0,0 +1,80 @@
+# Element types other than Float64. The internal convergence tolerances are
+# multiples of `eps(T)`, so both a narrower and a wider type must behave: a
+# Float64-scaled literal is unreachable in Float32 (every descent would run to
+# `maxiter`) and stops far short of what BigFloat can resolve.
+
+@testset "element types" begin
+
+    @testset "Float32 flows through the family" begin
+        A = Float32[4 1.5; 1.5 1]
+        for a in (symcover(A), soft_symcover(A), soft_symcover(AbsLog{1}(), A),
+                  soft_symcover(AbsLinear{1}(), A), symcover_min(AbsLog{2}(), A),
+                  soft_symcover_min(AbsLog{2}(), A))
+            @test eltype(a) === Float32
+            @test all(isfinite, a)
+        end
+        @test iscover(symcover(A), A; rtol=8eps(Float32))
+        @test iscover(symcover_min(AbsLog{2}(), A), A; rtol=8eps(Float32))
+
+        B = Float32[4 1.5 0.5; 1.5 1 2]
+        a, b = cover(B)
+        @test eltype(a) === eltype(b) === Float32
+        @test iscover(a, b, B; rtol=8eps(Float32))
+    end
+
+    @testset "BigFloat flows through the family" begin
+        A = BigFloat[4 1.5; 1.5 1]
+        for a in (symcover(A), soft_symcover(A), symcover_min(AbsLog{2}(), A),
+                  soft_symcover_min(AbsLog{2}(), A))
+            @test eltype(a) === BigFloat
+            @test all(isfinite, a)
+        end
+        @test iscover(symcover(A), A; rtol=8eps(BigFloat))
+    end
+
+    @testset "row and column scales may differ in element type" begin
+        # The public entry points take plain `AbstractVector`s, so a mismatch must
+        # not surface as a MethodError from an unexported internal. Each vector
+        # keeps its own element type; the arithmetic promotes.
+        A = [4.0 1.5 0.5; 1.5 1.0 2.0]
+        a0, b0 = Float64[3.0, 3.0], Float32[3.0, 3.0, 3.0]
+
+        a, b = cover!(copy(a0), copy(b0), A)
+        @test eltype(a) === Float64
+        @test eltype(b) === Float32
+        @test iscover(a, b, A; rtol=1e-6)
+
+        for ϕ in (AbsLog{1}(), AbsLinear{1}(), AbsLinear{2}())
+            x, y = soft_cover!(ϕ, copy(a0), copy(b0), A)
+            @test eltype(x) === Float64
+            @test eltype(y) === Float32
+            @test all(isfinite, x) && all(isfinite, y)
+        end
+
+        # The internals the public path routes through, reached directly.
+        @test MatrixCovers.tighten_cover!(copy(a0), copy(b0), A) isa Tuple
+        @test MatrixCovers.tighten_cover!(copy(b0), copy(a0), A') isa Tuple
+        @test MatrixCovers.unconstrained_min!(AbsLog{2}(), copy(a0), copy(b0), A) isa Tuple
+    end
+
+    @testset "tolerances follow the element type" begin
+        # Each convergence threshold must sit above the type's resolution, or the
+        # test it guards can never fire.
+        for T in (Float32, Float64, BigFloat)
+            @test MatrixCovers._multistart_switchtol(T) > eps(T)
+        end
+
+        # And must sit below it for a wider type, so the extra precision is used.
+        # Running the ALS kernel from one start under the eltype-scaled tolerance
+        # and under the Float64-scaled literal it replaced separates the two.
+        A = BigFloat[4 1.5 0.3; 1.5 1 0.7; 0.3 0.7 2.0]
+        start = initialize_symcover(A; feasible=:none)
+        u1, v1 = copy(start), copy(start)
+        u2, v2 = copy(start), copy(start)
+        MatrixCovers._mscm_als!(u1, v1, A, 500)
+        MatrixCovers._mscm_als!(u2, v2, A, 500; tol=1e-14)
+        @test cover_objective(AbsLinear{2}(), u1, v1, A) <
+              cover_objective(AbsLinear{2}(), u2, v2, A)
+    end
+
+end
diff --git a/test/extensions.jl b/test/extensions.jl
index aaf5235..a0349f5 100644
--- a/test/extensions.jl
+++ b/test/extensions.jl
@@ -49,6 +49,18 @@
     A_rank1 = [2.0 1.0 4.0; 1.0 0.5 2.0; 4.0 2.0 8.0]
     a_soft = soft_symcover_min(AbsLog{2}(), A_rank1)
     @test cover_objective(AbsLog{2}(), a_soft, A_rank1) < 1e-8
+
+    # A solve that does not reach an optimum is an error, not a cover: the point
+    # such a model holds is the base of a ray rather than a minimizer. The guard is
+    # exercised directly on the status, since the symmetry precondition rejects the
+    # asymmetric input that is what left this LP unbounded.
+    @test MatrixCovers.check_solved(JuMP.OPTIMAL, "HiGHS", "symcover_min") === nothing
+    @test MatrixCovers.check_solved(JuMP.LOCALLY_SOLVED, "Ipopt", "symcover_min!") === nothing
+    @test_throws "terminated with status" MatrixCovers.check_solved(JuMP.DUAL_INFEASIBLE, "HiGHS", "symcover_min")
+    @test_throws "symcover_min" MatrixCovers.check_solved(JuMP.DUAL_INFEASIBLE, "HiGHS", "symcover_min")
+    @test_throws "HiGHS" MatrixCovers.check_solved(JuMP.INFEASIBLE, "HiGHS", "symcover_min")
+    # A tolerance the caller did not ask for is not a solve.
+    @test_throws "terminated with status" MatrixCovers.check_solved(JuMP.ALMOST_LOCALLY_SOLVED, "Ipopt", "symcover_min!")
 end
 
 @testset "symcover_min and soft_symcover_min (JuMP/Ipopt, AbsLinear)" begin
@@ -140,16 +152,16 @@ end
     # The AbsLinear objectives are non-convex: on this matrix the :hardcover and
     # :geomean starts descend into genuinely different local minima, which is what
     # makes a menu of starts worth having.
-    Abasin = [6.609216272192496 1.032613546278995 55.276094662396076 0.5076927138328724;
-              1.032613546278995 3.1390835186570034 11.658167585612446 38.315826566607555;
-              55.276094662396076 11.658167585612446 0.001705708114264713 21.68951642774627;
-              0.5076927138328724 38.315826566607555 21.68951642774627 0.006443251375371587]
+    Abasin = [81.892035218799 1.06622031288736 29.4700945830419 0.0181293142917846;
+              1.06622031288736 0.243512973596586 38.0236584552296 0.0279078887878805;
+              29.4700945830419 38.0236584552296 8.96405068596511 26.5775238859338;
+              0.0181293142917846 0.0279078887878805 26.5775238859338 42.6650094474717]
     ϕ = AbsLinear{2}()
     a_hard = symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin)
     a_geo = symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:geomean), Abasin)
     @test iscover(a_hard, Abasin; rtol=1e-6)
     @test iscover(a_geo, Abasin; rtol=1e-6)
-    @test cover_objective(ϕ, a_geo, Abasin) < cover_objective(ϕ, a_hard, Abasin) - 1e-3
+    @test cover_objective(ϕ, a_geo, Abasin) < cover_objective(ϕ, a_hard, Abasin) - 0.5
 
     # Offset axes survive the Ipopt position mapping.
     Ao = OffsetArray(A, -1, -1)
@@ -300,10 +312,10 @@ end
     # The matrix on which the starts genuinely disagree: :geomean reaches a better
     # AbsLinear{2} minimum than :hardcover, so a driver that tried only the latter
     # would report the worse of the two.
-    Abasin = [6.609216272192496 1.032613546278995 55.276094662396076 0.5076927138328724;
-              1.032613546278995 3.1390835186570034 11.658167585612446 38.315826566607555;
-              55.276094662396076 11.658167585612446 0.001705708114264713 21.68951642774627;
-              0.5076927138328724 38.315826566607555 21.68951642774627 0.006443251375371587]
+    Abasin = [81.892035218799 1.06622031288736 29.4700945830419 0.0181293142917846;
+              1.06622031288736 0.243512973596586 38.0236584552296 0.0279078887878805;
+              29.4700945830419 38.0236584552296 8.96405068596511 26.5775238859338;
+              0.0181293142917846 0.0279078887878805 26.5775238859338 42.6650094474717]
     Aasym = [1.0 2.0 3.0; 40.0 5.0 0.6]
 
     for ϕ in (AbsLinear{1}(), AbsLinear{2}())
@@ -339,7 +351,7 @@ end
     @test symcover_min(ϕ, Abasin; strategies=(:hardcover,)) ≈
           symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin)
     @test cover_objective(ϕ, symcover_min(ϕ, Abasin), Abasin) <
-          cover_objective(ϕ, symcover_min(ϕ, Abasin; strategies=(:hardcover,)), Abasin) - 1e-3
+          cover_objective(ϕ, symcover_min(ϕ, Abasin; strategies=(:hardcover,)), Abasin) - 0.5
 
     # A matrix whose every row carries a single support entry admits no :leaveout start;
     # that strategy forfeits its slot rather than failing the solve.
@@ -349,7 +361,14 @@ end
     @test_throws "unknown strategy :banana" symcover_min(ϕ, Abasin; strategies=(:banana,))
     @test_throws "no strategy in (:leaveout,)" symcover_min(ϕ, Adiag; strategies=(:leaveout,))
     @test_throws "has no asymmetric formulation" cover_min(ϕ, Aasym; strategies=(:leaveout,))
+    # An empty menu is reported as such by all four drivers, rather than as the
+    # unrelated "no strategy yields a start" that a fall-through would produce.
     @test_throws "at least one starting cover" cover_min(ϕ, Aasym; strategies=())
+    @test_throws "at least one starting cover" symcover_min(ϕ, Abasin; strategies=())
+    @test_throws "at least one starting cover" soft_cover_min(ϕ, Aasym; strategies=())
+    @test_throws "at least one starting cover" soft_symcover_min(ϕ, Abasin; strategies=())
+    @test_throws "symcover_min:" symcover_min(ϕ, Abasin; strategies=())
+    @test_throws "soft_symcover_min:" soft_symcover_min(ϕ, Abasin; strategies=())
 end
 
 @testset "error hint gated on argument types" begin
@@ -454,3 +473,64 @@ end
     @test e5 isa MethodError
     @test !occursin("loading JuMP", sprint(showerror, e5))
 end
+
+# `HookOnlyMatrix` (defined in soft_covers.jl) throws from `getindex`, so a model
+# builder that materialized `A` in position space could not get this far.
+@testset "the solver entry points read through the support hook" begin
+    entries = [(1, 1, 2.0), (1, 3, 1.5), (2, 2, 3.0), (2, 4, 0.5), (3, 4, 4.0), (4, 4, 1.0)]
+    M = HookOnlyMatrix(entries, 4)
+    dense = zeros(4, 4)
+    for (i, j, v) in entries
+        dense[i, j] = dense[j, i] = v
+    end
+
+    @test symcover_min(AbsLog{2}(), M) ≈ symcover_min(AbsLog{2}(), dense) rtol=1e-8
+    @test MatrixCovers.symcover_min_jump(AbsLog{2}(), M) ≈
+          MatrixCovers.symcover_min_jump(AbsLog{2}(), dense) rtol=1e-6
+    @test symcover_min(AbsLog{1}(), M) ≈ symcover_min(AbsLog{1}(), dense) rtol=1e-6
+    @test iscover(symcover_min(AbsLog{1}(), M), dense; rtol=1e-8)
+
+    for (fh, fd) in ((cover_min(AbsLog{1}(), M), cover_min(AbsLog{1}(), dense)),
+                     (MatrixCovers.cover_min_jump(AbsLog{2}(), M),
+                      MatrixCovers.cover_min_jump(AbsLog{2}(), dense)))
+        @test fh[1] ≈ fd[1] rtol=1e-6
+        @test fh[2] ≈ fd[2] rtol=1e-6
+    end
+
+    # Ipopt: the AbsLinear kernels, whose starts the caller supplies.
+    for p in (1, 2)
+        ϕ = AbsLinear{p}()
+        sh = symcover_min!(ϕ, symcover(ϕ, M), M)
+        sd = symcover_min!(ϕ, symcover(ϕ, dense), dense)
+        @test sh ≈ sd rtol=1e-5
+
+        ah, bh = cover_min!(ϕ, cover(ϕ, M)..., M)
+        ad, bd = cover_min!(ϕ, cover(ϕ, dense)..., dense)
+        @test ah ≈ ad rtol=1e-5
+        @test bh ≈ bd rtol=1e-5
+
+        @test soft_symcover_min!(ϕ, soft_symcover(ϕ, M), M) ≈
+              soft_symcover_min!(ϕ, soft_symcover(ϕ, dense), dense) rtol=1e-5
+
+        sah, sbh = soft_cover_min!(ϕ, soft_cover(ϕ, M)..., M)
+        sad, sbd = soft_cover_min!(ϕ, soft_cover(ϕ, dense)..., dense)
+        @test sah ≈ sad rtol=1e-5
+        @test sbh ≈ sbd rtol=1e-5
+    end
+end
+
+# The `sym` AbsLinear solvers minimize the full-grid objective — each off-diagonal
+# pair twice, each diagonal entry once — which is what `cover_objective` reports.
+# This matrix discriminates between that and weighting each unordered pair once:
+# the latter convention puts the optimum near [2.0, 3.177, 1.574] instead. Most
+# matrices do not discriminate, because the binding constraints have zero residual
+# at the optimum and a zero residual is weight-independent.
+@testset "the Ipopt sym objective uses the full-grid weighting" begin
+    A = Float64[4 1 0; 1 1 5; 0 5 2]
+    a0 = fill(sqrt(maximum(A)), 3)
+    a = symcover_min!(AbsLinear{2}(), copy(a0), A)
+    @test a ≈ [2.0, 1.0, 5.0] rtol=1e-6
+    @test iscover(a, A; rtol=1e-6)
+    # The objective reported for the result is the one that was minimized.
+    @test cover_objective(AbsLinear{2}(), a, A) ≈ cover_objective(AbsLinear{2}(), [2.0, 1.0, 5.0], A) rtol=1e-6
+end
diff --git a/test/heuristic_covers.jl b/test/heuristic_covers.jl
index 2808c89..e3bec33 100644
--- a/test/heuristic_covers.jl
+++ b/test/heuristic_covers.jl
@@ -45,6 +45,18 @@ end
             @test symcover(ϕ, A) == a
         end
     end
+
+    # Ignored is not the same as unchecked: the slot takes a penalty, so a wrong
+    # first argument fails rather than being silently dropped.
+    A = [4.0 1.5; 1.5 1.0]
+    a = symcover(A)
+    b = copy(a)
+    for bad in (42, "nope", :whatever)
+        @test_throws MethodError symcover(bad, A)
+        @test_throws MethodError symcover!(bad, copy(a), A)
+        @test_throws MethodError cover(bad, A)
+        @test_throws MethodError cover!(bad, copy(a), copy(b), A)
+    end
 end
 
 @testset "symcover with unequal row degrees" begin
@@ -207,13 +219,13 @@ end
     @test median(gaps) < 0.0195 * 1.5
 end
 
-@testset "boost feasibility: lower-dominant bands and extreme dynamic range" begin
-    # Lower-dominant bands: the symmetric-contract value is the max over both
-    # triangles, so the subdiagonal must be covered even when it dominates.
-    a = symcover(AbsLog{2}(), Bidiagonal([0.01, 0.01, 0.01], [100.0, 100.0], 'L'); maxiter=0)
+@testset "boost feasibility: dominant off-diagonal bands and extreme dynamic range" begin
+    # An off-diagonal band far larger than the diagonal: the boost must raise both
+    # endpoints of each band entry, not just the diagonal it started from.
+    a = symcover(AbsLog{2}(), SymTridiagonal([0.01, 0.01, 0.01], [100.0, 100.0]); maxiter=0)
     @test a[1] * a[2] >= 100 * (1 - 8eps())
     @test a[2] * a[3] >= 100 * (1 - 8eps())
-    T40 = Tridiagonal(fill(50.0, 39), fill(0.01, 40), fill(0.02, 39))
+    T40 = Tridiagonal(fill(50.0, 39), fill(0.01, 40), fill(50.0, 39))
     a = symcover(AbsLog{2}(), T40; maxiter=0)
     M = Matrix(T40)
     @test iscover(a, M; rtol=8eps())
diff --git a/test/invariants.jl b/test/invariants.jl
index bf47f35..8d6e1a9 100644
--- a/test/invariants.jl
+++ b/test/invariants.jl
@@ -96,4 +96,53 @@ const GEN_NOTIONS = (
             strategy in (:hardcover, :geomean), feasible in (:inflate, :boost, :none)
         @test isbalanced(initialize_cover(Agen; strategy, feasible)..., Agen)
     end
+
+    # The sym objective is summed over the full grid: each off-diagonal pair counts
+    # twice, each diagonal entry once. `cover_objective` is the reference, and every
+    # sym solver must minimize that same weighting even though its constraints live on
+    # the `i <= j` triangle. The references below impose the constraints on the full
+    # grid explicitly, so agreeing with them pins both halves of the convention: that
+    # the triangle is the equivalent constraint set, and that the objective is not
+    # halved along with it.
+    @testset "sym solvers minimize the full-grid objective" begin
+        function ref_min(pow, A)
+            n = size(A, 1)
+            L = log.(abs.(A))
+            supp = [(i, j) for i in 1:n, j in 1:n if !iszero(A[i, j])]
+            model = JuMP.Model(HiGHS.Optimizer)
+            JuMP.set_silent(model)
+            JuMP.@variable(model, α[1:n])
+            r(i, j) = α[i] + α[j] - L[i, j]
+            JuMP.@objective(model, Min, sum(r(i, j)^pow for (i, j) in supp))
+            for (i, j) in supp
+                JuMP.@constraint(model, r(i, j) >= 0)
+            end
+            JuMP.optimize!(model)
+            return JuMP.objective_value(model)
+        end
+
+        # The two conventions share a minimizer whenever the diagonal residuals
+        # vanish at the optimum, which is the common case — so a discriminating
+        # matrix is committed rather than left to the random draws. Weighting this
+        # one's off-diagonal pairs once instead of twice moves the optimum.
+        Mdisc = [2.4 1.2 1.2; 1.2 1.4 2.4; 1.2 2.4 0.6]
+
+        rng = StableRNG(20)
+        mats = Any[Asym, Azsym, abs.(Hc), Mdisc]
+        for _ in 1:5
+            n = rand(rng, 3:7)
+            B = randn(rng, n, n) .* (rand(rng, n, n) .< 0.5)
+            M = B + transpose(B)
+            push!(mats, M)
+        end
+        for M in mats
+            # AbsLog{2}: native solver and the HiGHS reference model.
+            o2 = ref_min(2, M)
+            @test cover_objective(AbsLog{2}(), symcover_min(AbsLog{2}(), M), M) ≈ o2 rtol=1e-5
+            @test cover_objective(AbsLog{2}(), MatrixCovers.symcover_min_jump(AbsLog{2}(), M), M) ≈ o2 rtol=1e-5
+            # AbsLog{1}: the HiGHS LP, whose objective is assembled from nonzero counts
+            # rather than from the residuals directly.
+            @test cover_objective(AbsLog{1}(), symcover_min(AbsLog{1}(), M), M) ≈ ref_min(1, M) rtol=1e-5
+        end
+    end
 end
diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl
index d254b73..31e837d 100644
--- a/test/minimal_covers.jl
+++ b/test/minimal_covers.jl
@@ -317,3 +317,24 @@ end
     end
     @test median(gen_ratios) < 1.02
 end
+
+# The `:lsqr` path is the intended solve when nnz ≪ n²; its working set must
+# therefore be sized by the support, not by `length(A)`. Doubling `n` at fixed
+# nnz-per-row doubles the support, so allocation must roughly double too — a
+# working set carrying any n×n array would quadruple instead.
+@testset ":lsqr allocates in proportion to the support" begin
+    function lsqr_alloc(n)
+        rng = StableRNG(4)
+        B = sprandn(rng, n, n, 3.5 / n)
+        A = Symmetric(B + transpose(B))
+        MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr)   # compile before measuring
+        return @allocated MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr)
+    end
+    small, large = lsqr_alloc(200), lsqr_alloc(800)
+    # Measured ratio is ≈6: the support grows 4x, and the Newton continuation takes
+    # more LSQR iterations at the larger size. Quadratic growth would put it near 16.
+    # The bound discriminates sharply despite the gap, because one n×n Float64 array
+    # at n = 800 is 4.9 MB against a whole measured working set of ≈3.9 MB — a single
+    # reintroduced one lands the ratio above 13.
+    @test large < 10 * small
+end
diff --git a/test/penalties.jl b/test/penalties.jl
index 0f80cfc..ac479e3 100644
--- a/test/penalties.jl
+++ b/test/penalties.jl
@@ -23,6 +23,64 @@
     @test cover_objective(AbsLinear{2}(), a0, A0) ≈ 0.0 + 1.0 + 1.0 + 1/16   # (0,1), (1,0) off-diag zeros
 end
 
+# A matrix readable only through the traversal hook: `getindex` throws, so any
+# full-grid scan fails outright.
+struct SupportOnlyMatrix{T} <: AbstractMatrix{T}
+    entries::Vector{Tuple{Int,Int,T}}
+    sz::Tuple{Int,Int}
+end
+Base.size(M::SupportOnlyMatrix) = M.sz
+Base.getindex(::SupportOnlyMatrix, ::Int, ::Int) =
+    error("SupportOnlyMatrix must be read through foreach_support")
+function MatrixCovers.foreach_support(f, M::SupportOnlyMatrix)
+    for (i, j, v) in M.entries
+        iszero(v) || f(i, j, abs(v))
+    end
+    return nothing
+end
+
+@testset "cover_objective reads through the support hook" begin
+    entries = [(1, 2, 3.0), (2, 1, -1.5), (3, 3, 4.0)]
+    M = SupportOnlyMatrix(entries, (3, 4))
+    dense = zeros(3, 4)
+    for (i, j, v) in entries
+        dense[i, j] = v
+    end
+    a, b = [2.0, 1.0, 0.5], [1.5, 0.5, 3.0, 1.0]
+    for ϕ in PENALTIES
+        # Agreement with the dense reference pins the zero-entry accounting: the
+        # entries the hook skips still carry `ϕ(0)`, which `AbsLinear` makes nonzero.
+        @test cover_objective(ϕ, a, b, M) ≈ cover_objective(ϕ, a, b, dense)
+    end
+
+    # A zero scale under a nonzero entry is uncovered, and under a zero entry is not.
+    az = [0.0, 1.0, 0.5]
+    @test isinf(cover_objective(AbsLog{2}(), az, b, M))
+    # Column 4 carries no support, so zeroing its scale constrains nothing.
+    @test isfinite(cover_objective(AbsLog{2}(), a, [1.5, 0.5, 3.0, 0.0], M))
+
+    @test cover_objective(AbsLog{2}(), b, a, M') ≈ cover_objective(AbsLog{2}(), a, b, dense)
+
+    # A penalty infinite at `r = 0` is legal, and a matrix with no zero entries must
+    # not pick up its value: the zero-entry term is skipped, not multiplied by zero.
+    infatzero(r) = iszero(r) ? Inf : abs(log(r))
+    @test cover_objective(infatzero, [1.0, 1.0], [1.0, 1.0], [1.0 2.0; 3.0 4.0]) ≈
+          sum(abs(log(v)) for v in (1.0, 2.0, 3.0, 4.0))
+    @test isinf(cover_objective(infatzero, [1.0, 1.0], [1.0, 1.0], [1.0 0.0; 3.0 4.0]))
+end
+
+@testset "cover_objective: index checking" begin
+    A = [4.0 1.5; 1.5 1.0]
+    @test_throws "indices of `a` must match row-indexing of `A`" cover_objective(AbsLog{2}(), [1.0], [1.0, 1.0], A)
+    @test_throws "indices of `b` must match column-indexing of `A`" cover_objective(AbsLog{2}(), [1.0, 1.0], [1.0], A)
+    # Offset axes are honored, not merely tolerated: the score is unchanged.
+    Ao = OffsetArray(A, -1:0, 2:3)
+    ao, bo = OffsetArray([2.0, 1.0], -1:0), OffsetArray([1.0, 3.0], 2:3)
+    for ϕ in PENALTIES
+        @test cover_objective(ϕ, ao, bo, Ao) ≈ cover_objective(ϕ, [2.0, 1.0], [1.0, 3.0], A)
+    end
+end
+
 @testset "cover_objective: complex input" begin
     # The objective depends only on entry magnitudes, so complex A and abs.(A)
     # give identical results, and the accumulator stays real.
diff --git a/test/runtests.jl b/test/runtests.jl
index 1e43759..c6d33d2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,6 +24,7 @@ include("helpers.jl")               # isbalanced, covaries, PENALTIES
     include("initializers.jl")      # initialize_symcover/initialize_cover strategies
     include("minimal_covers.jl")    # the *_min family (native solvers)
     include("storage_types.jl")     # sparse/structured/wrapped storage vs dense reference
+    include("element_types.jl")     # Float32/BigFloat: eltype-scaled internal tolerances
     include("extensions.jl")        # JuMP/HiGHS and Ipopt solvers, missing-extension hints
     include("unitful.jl")           # dimensional covers via the Unitful extension
     include("invariants.jl")        # shared conventions checked across every notion
@@ -42,12 +43,16 @@ include("helpers.jl")               # isbalanced, covaries, PENALTIES
                      :_prepare_cover_start!, :_prepare_symcover_start!,
                      :_prepare_soft_cover_start!, :_prepare_soft_symcover_start!,
                      :foreach_support, :foreach_support_sym,
-                     :cover_min_jump, :symcover_min_jump)
+                     :cover_min_jump, :symcover_min_jump, :check_solved,
+                     :require_abs_symmetric,
+                     :_edge_list, :_sym_edge_list, :_degrees)
         # Non-public names owned by other packages, each with no public equivalent:
-        # `FreeUnits`/`Unit` are Unitful's unit representation, `Optimizer` is the
-        # solver handle JuMP's own documented `Model(HiGHS.Optimizer)` entry point
-        # names, and `register_error_hint` is Base-internal.
-        foreign = (:FreeUnits, :Unit, :Optimizer, :Experimental, :register_error_hint)
+        # `FreeUnits`/`Unit` are Unitful's unit representation and `Units` their
+        # abstract supertype, needed to reject the unit types this package cannot
+        # read; `Optimizer` is the solver handle JuMP's own documented
+        # `Model(HiGHS.Optimizer)` entry point names; and `register_error_hint` is
+        # Base-internal.
+        foreign = (:FreeUnits, :Unit, :Units, :Optimizer, :Experimental, :register_error_hint)
         test_explicit_imports(
             MatrixCovers;
             all_explicit_imports_are_public = VERSION >= v"1.11" ?
diff --git a/test/soft_covers.jl b/test/soft_covers.jl
index 5916d4c..03413d1 100644
--- a/test/soft_covers.jl
+++ b/test/soft_covers.jl
@@ -526,3 +526,92 @@ end
         @test collect(ao) ≈ aref rtol=1e-6
     end
 end
+
+# A matrix readable only through the support hook: any full-grid scan hits the
+# throwing `getindex`. Both traversals are defined so the sym and asym kernels
+# can be driven from the same fixture.
+struct HookOnlyMatrix{T} <: AbstractMatrix{T}
+    entries::Vector{Tuple{Int,Int,T}}   # `i <= j` only; the transpose is implied
+    n::Int
+end
+Base.size(M::HookOnlyMatrix) = (M.n, M.n)
+Base.getindex(::HookOnlyMatrix, ::Int, ::Int) =
+    error("HookOnlyMatrix must be read through foreach_support")
+function MatrixCovers.foreach_support_sym(f, M::HookOnlyMatrix)
+    for (i, j, v) in M.entries
+        iszero(v) || f(i, j, abs(v))
+    end
+    return nothing
+end
+function MatrixCovers.foreach_support(f, M::HookOnlyMatrix)
+    for (i, j, v) in M.entries
+        iszero(v) && continue
+        f(i, j, abs(v))
+        i == j || f(j, i, abs(v))
+    end
+    return nothing
+end
+# Storing one member of each pair makes `abs`-symmetry structural, so the
+# precondition check is a no-op — as it must be, since it cannot index `M`.
+MatrixCovers.require_abs_symmetric(::HookOnlyMatrix, fname) = nothing
+
+@testset "the soft-cover kernels read through the support hook" begin
+    entries = [(1, 1, 2.0), (1, 3, 1.5), (2, 2, 3.0), (2, 4, 0.5), (3, 4, 4.0), (4, 4, 1.0)]
+    M = HookOnlyMatrix(entries, 4)
+    dense = zeros(4, 4)
+    for (i, j, v) in entries
+        dense[i, j] = dense[j, i] = v
+    end
+    start() = [1.3, 0.7, 2.1, 0.9]
+
+    # Reaching a result at all proves the kernel never indexed `M`; matching the
+    # dense run proves the gathered support is the same set of entries carrying the
+    # same full-grid multiplicity.
+    for kernel! in (MatrixCovers._abslog1_iter!, MatrixCovers._abslinear1_iter!,
+                    MatrixCovers._abslinear2_iter!)
+        @test kernel!(start(), M, 50) ≈ kernel!(start(), dense, 50) rtol=1e-10
+    end
+    for kernel! in (MatrixCovers._abslog1_iter_asym!, MatrixCovers._abslinear1_iter_asym!,
+                    MatrixCovers._mscm_als!)
+        ah, bh = kernel!(start(), start(), M, 50)
+        ad, bd = kernel!(start(), start(), dense, 50)
+        @test ah ≈ ad rtol=1e-10
+        @test bh ≈ bd rtol=1e-10
+    end
+end
+
+# The kernels above are only half the path: the public entry points reach them
+# through the initializers, so a full-grid scan in either one would surface here.
+@testset "the cover entry points read through the support hook" begin
+    entries = [(1, 1, 2.0), (1, 3, 1.5), (2, 2, 3.0), (2, 4, 0.5), (3, 4, 4.0), (4, 4, 1.0)]
+    M = HookOnlyMatrix(entries, 4)
+    dense = zeros(4, 4)
+    for (i, j, v) in entries
+        dense[i, j] = dense[j, i] = v
+    end
+
+    # `:diagfeasible` runs `init_feasible_diag!` and `boost_feasible_seq!`;
+    # `:leaveout` runs `_leaveout_logmean_init!`.
+    for strategy in (:geomean, :hardcover, :leaveout, :diagfeasible)
+        for feasible in (:none, :boost, :inflate)
+            @test initialize_symcover(M; strategy, feasible) ≈
+                  initialize_symcover(dense; strategy, feasible) rtol=1e-10
+        end
+    end
+
+    for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}(), AbsLog{2}())
+        @test symcover(ϕ, M) ≈ symcover(ϕ, dense) rtol=1e-10
+        @test soft_symcover(ϕ, M) ≈ soft_symcover(ϕ, dense) rtol=1e-6
+        @test iscover(symcover(ϕ, M), dense; rtol=8eps())
+        @test cover_objective(ϕ, symcover(ϕ, M), M) ≈
+              cover_objective(ϕ, symcover(ϕ, dense), dense) rtol=1e-10
+    end
+
+    for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}(), AbsLog{2}())
+        ah, bh = cover(ϕ, M)
+        ad, bd = cover(ϕ, dense)
+        @test ah ≈ ad rtol=1e-10
+        @test bh ≈ bd rtol=1e-10
+        @test iscover(ah, bh, dense; rtol=8eps())
+    end
+end
diff --git a/test/storage_types.jl b/test/storage_types.jl
index b6509d6..05975a1 100644
--- a/test/storage_types.jl
+++ b/test/storage_types.jl
@@ -1,6 +1,28 @@
 # Storage-type equivalence: sparse, structured, and wrapped inputs must match the
 # dense reference.
 
+@testset "complex Hermitian sparse" begin
+    # Only `abs` of a stored value is read, so a complex `Hermitian` gets the
+    # sparse traversal rather than falling back to the generic O(n^2) getindex.
+    P = sparse([1, 2, 1, 3], [1, 2, 2, 3], ComplexF64[4.0, 1.0, 0.5+0.5im, 2.0], 3, 3)
+    for uplo in (:U, :L)
+        H = Hermitian(P, uplo)
+        @test which(MatrixCovers.foreach_support_sym, Tuple{Function,typeof(H)}).file ==
+              which(MatrixCovers.foreach_support_sym,
+                    Tuple{Function,Symmetric{Float64,SparseMatrixCSC{Float64,Int}}}).file
+
+        got = Tuple{Int,Int,Float64}[]
+        MatrixCovers.foreach_support_sym((i, j, v) -> push!(got, (i, j, v)), H)
+        ref = Tuple{Int,Int,Float64}[]
+        MatrixCovers.foreach_support_sym((i, j, v) -> push!(ref, (i, j, v)), Matrix(H))
+        @test sort(got) == sort(ref)
+
+        a = symcover_min(AbsLog{2}(), H)
+        @test iscover(a, Matrix(H); rtol=1e-8)
+        @test symcover(AbsLog{2}(), H) ≈ symcover(AbsLog{2}(), Matrix(H)) rtol = 1e-12
+    end
+end
+
 @testset "SparseMatrixCSC" begin
     for Adense in ([2.0 1.0; 1.0 3.0], [1.0 -0.2; -0.2 0.0], [0.0 12.0 9.0; 12.0 7.0 12.0; 9.0 12.0 0.0],
                    [100.0 1.0; 1.0 0.01])
@@ -243,3 +265,26 @@ end
               cover_objective(AbsLog{2}(), cover_min(AbsLog{2}(), M)..., M) rtol = 1e-7
     end
 end
+
+# `Symmetric`/`Hermitian` over a sparse parent store one triangle, so the
+# asymmetric traversal must reconstitute the other. Without its own method these
+# fall back to the generic full-grid `getindex` scan, which is correct but defeats
+# the sparse specialization the wrapper exists to enable.
+@testset "asymmetric traversal of wrapped sparse storage" begin
+    P = sparse([1, 2, 1, 3], [1, 2, 3, 3], [2.0, 3.0, 1.5, 4.0], 3, 3)
+    for W in (Symmetric(P, :U), Symmetric(sparse(transpose(P)), :L),
+              Hermitian(complex(P), :U))
+        emitted = Tuple{Int,Int,Float64}[]
+        foreach_support(W) do i, j, v
+            push!(emitted, (i, j, v))
+        end
+        reference = Tuple{Int,Int,Float64}[]
+        foreach_support(Matrix(W)) do i, j, v
+            push!(reference, (i, j, v))
+        end
+        @test sort(emitted) == sort(reference)
+        # Each pair exactly once per orientation: a doubled emission would silently
+        # double that entry's weight in every objective.
+        @test length(emitted) == length(unique(e -> (e[1], e[2]), emitted))
+    end
+end
diff --git a/test/support.jl b/test/support.jl
index 34b4db2..2a66ad0 100644
--- a/test/support.jl
+++ b/test/support.jl
@@ -97,3 +97,103 @@
         end
     end
 end
+
+@testset "the symmetry precondition" begin
+    # A genuinely asymmetric matrix is refused rather than covered as some
+    # unspecified symmetrization of itself.
+    for A in ([0.0 0.0; 1.0 0.0], [1.0 2.0; 3.0 1.0], [4.0 1.0; -1.5 4.0])
+        @test_throws "abs.(A)` to be symmetric" symcover(A)
+        @test_throws "abs.(A)` to be symmetric" soft_symcover(A)
+        @test_throws "abs.(A)` to be symmetric" symcover_min(AbsLog{2}(), A)
+        @test_throws "abs.(A)` to be symmetric" initialize_symcover(A)
+        @test_throws "abs.(A)` to be symmetric" sparse(A) |> symcover
+        # Every solver family, including the ones behind the extensions.
+        @test_throws "abs.(A)` to be symmetric" symcover_min(AbsLog{1}(), A)
+        @test_throws "abs.(A)` to be symmetric" symcover_min(AbsLinear{2}(), A)
+        @test_throws "abs.(A)` to be symmetric" soft_symcover_min(AbsLinear{2}(), A)
+    end
+    # The message names the offending pair and the entry point.
+    @test_throws "symcover!" symcover([0.0 0.0; 1.0 0.0])
+    @test_throws "abs(A[2,1])" symcover([0.0 0.0; 1.0 0.0])
+
+    # Only the magnitudes must agree, so a sign flip is fine and a complex
+    # Hermitian — where A[i,j] == conj(A[j,i]) — is a legitimate input.
+    @test iscover(symcover([4.0 1.0; -1.0 4.0]), [4.0 1.0; -1.0 4.0]; rtol=8eps())
+    P = sparse([1, 2, 1], [1, 2, 2], ComplexF64[4.0, 2.0, 0.5+0.5im], 2, 2)
+    @test symcover_min(AbsLog{2}(), Hermitian(P, :U)) isa AbstractVector
+
+    # Exact symmetry is not required: a symmetric matrix that has been through
+    # floating-point arithmetic lands a ULP or so off, and must still be accepted.
+    rng = StableRNG(7)
+    B = randn(rng, 12, 12); S = B + B'
+    d = exp.(randn(rng, 12))
+    Sd = (d .* S) .* d'          # symmetric in exact arithmetic only
+    @test !issymmetric(Sd)
+    @test symcover(Sd) isa AbstractVector
+
+    # Wrapper types are exempt because their storage makes the property structural.
+    M = [4.0 1.0; 2.0 4.0]       # asymmetric as written
+    @test symcover(Symmetric(M, :U)) isa AbstractVector
+    @test symcover(Symmetric(M, :L)) isa AbstractVector
+    @test symcover(Diagonal([1.0, 2.0])) isa AbstractVector
+
+    # Banded storage earns no exemption. A `Bidiagonal` reads one of its
+    # off-diagonals as a structural zero, so any nonzero band makes it asymmetric;
+    # a `Tridiagonal` stores both bands and qualifies only when they agree.
+    for A in (Bidiagonal([3.0, 2.0, 1.0], [6.0, 0.5], :U),
+              Bidiagonal([3.0, 2.0, 1.0], [6.0, 0.5], :L),
+              Tridiagonal([1.0, 0.5], [3.0, 2.0, 1.0], [4.0, 0.5]))
+        @test_throws "abs.(A)` to be symmetric" symcover(A)
+        @test_throws "abs.(A)` to be symmetric" symcover_min(AbsLog{2}(), A)
+        # The asymmetric family covers them as stored.
+        @test iscover(cover(A)..., Matrix(A); rtol=8eps())
+    end
+    # Bands that do agree are ordinary symmetric input.
+    @test iscover(symcover(Tridiagonal([2.0, 0.5], [3.0, 2.0, 1.0], [2.0, 0.5])),
+                  Matrix(Tridiagonal([2.0, 0.5], [3.0, 2.0, 1.0], [2.0, 0.5])); rtol=8eps())
+    @test symcover(Bidiagonal([3.0, 2.0, 1.0], [0.0, 0.0], :U)) isa AbstractVector
+end
+
+# The kernels that gather the support into per-group neighbor lists read those
+# lists in place of the matrix, so the gather must reproduce the matrix exactly —
+# including multiplicity. `_sym_support` in particular enters each off-diagonal
+# pair in both orientations, which is what gives a kernel accumulating over its
+# groups the full-grid weighting of `cover_objective` (each off-diagonal pair
+# twice, each diagonal entry once) rather than a halved one.
+@testset "grouped support reproduces the matrix" begin
+    function regroup(S, groups_are_rows::Bool, sz)
+        R = zeros(Float64, sz)
+        for g in S.ax, s in MatrixCovers._slots(S, g)
+            i, j = groups_are_rows ? (g, S.idx[s]) : (S.idx[s], g)
+            R[i, j] = S.val[s]
+        end
+        return R
+    end
+
+    mats = Any[[2.0 1.0 0.0; 1.0 0.0 3.0; 0.0 3.0 4.0],       # zeros on and off the diagonal
+               [0.0 0.0; 0.0 0.0],                             # empty support
+               Diagonal([2.0, 0.0, 5.0]),
+               SymTridiagonal([4.0, 3.0, 1.0], [2.0, 0.5])]
+    rng = StableRNG(7)
+    for _ in 1:4
+        n = rand(rng, 2:6)
+        B = randn(rng, n, n) .* (rand(rng, n, n) .< 0.5)
+        push!(mats, B + transpose(B))
+    end
+    for M in mats
+        S = MatrixCovers._sym_support(M, Float64)
+        @test regroup(S, true, size(M)) == abs.(Matrix(M))
+    end
+
+    # The asymmetric gathers need no multiplicity rule: one entry, one slot.
+    for M in Any[[1.0 0.0 2.0; 0.0 0.0 0.0; 3.0 4.0 0.0], randn(StableRNG(9), 4, 6)]
+        @test regroup(MatrixCovers._row_support(M, Float64), true, size(M)) == abs.(M)
+        @test regroup(MatrixCovers._col_support(M, Float64), false, size(M)) == abs.(M)
+    end
+
+    # Offset axes: groups are addressed by the matrix's own indices.
+    O = OffsetArray([2.0 1.0; 1.0 3.0], -1:0, -1:0)
+    SO = MatrixCovers._sym_support(O, Float64)
+    @test SO.ax == -1:0
+    @test sort([(SO.idx[s], SO.val[s]) for s in MatrixCovers._slots(SO, -1)]) == [(-1, 2.0), (0, 1.0)]
+end
diff --git a/test/unitful.jl b/test/unitful.jl
index 6e6907f..f9a1cfd 100644
--- a/test/unitful.jl
+++ b/test/unitful.jl
@@ -271,4 +271,65 @@
         @test_throws Unitful.DimensionError symcover_min!(AbsLog{2}(), a3, A)
     end
 
+    @testset "unit types the cover cannot use" begin
+        # An affine unit measures from a shifted origin, so no product a[i]*b[j]
+        # scales it: refuse it rather than silently reducing it to its atom.
+        Aaff = [1.0 2.0; 2.0 3.0]u"°C"
+        @test_throws "affine units are not supported" symcover(Aaff)
+        @test_throws "°C" symcover(Aaff)
+        @test_throws "affine units are not supported" cover(Aaff)
+
+        # ContextUnits/FixedUnits carry a conversion context this code does not
+        # read; they are named, not left to fail inside an unexported internal.
+        for U in (Unitful.ContextUnits(u"m", u"m"), Unitful.FixedUnits(u"m"))
+            Actx = [1.0 2.0; 2.0 3.0] .* U^2
+            @test_throws "unsupported unit type" symcover(Actx)
+            @test_throws "FreeUnits" symcover(Actx)
+        end
+
+        # The ordinary case is untouched.
+        @test unit.(symcover([4.0 1.5; 1.5 1.0]u"m^2")) == [u"m", u"m"]
+    end
+
+    @testset "scoring a dimensional cover" begin
+        # `unit(A[i,j]) == unit(a[i])*unit(b[j])` makes every ratio the objective
+        # forms dimensionless, so the score is a plain number and is identical to
+        # the score of the unit-stripped problem.
+        ϕs = (AbsLog{1}(), AbsLog{2}(), AbsLinear{1}(), AbsLinear{2}())
+
+        a = symcover(A)
+        for ϕ in ϕs
+            s = cover_objective(ϕ, a, A)
+            @test s isa Float64
+            @test s == cover_objective(ϕ, ustrip.(a), ustrip.(A))
+        end
+
+        # Heterogeneous units across coordinates, not just a single uniform unit.
+        H = Quantity[4.0u"m^-2"        1.0u"m^-1*mm^-1";
+                     1.0u"m^-1*mm^-1" 4.0u"mm^-2"]
+        h = symcover(H)
+        for ϕ in ϕs
+            @test cover_objective(ϕ, h, H) == cover_objective(ϕ, ustrip.(h), ustrip.(H))
+        end
+
+        # The asymmetric form, whose row and column scales carry different units.
+        ru, cu = (u"m", u"s"), (u"kg", u"K", u"m/s")
+        B = [1.0 / (r * c) for r in ru, c in cu] .* [1e3 1e-2 5.0; 2.0 1e4 1e-1]
+        ab, bb = cover(B)
+        for ϕ in ϕs
+            s = cover_objective(ϕ, ab, bb, B)
+            @test s isa Float64
+            @test s == cover_objective(ϕ, ustrip.(ab), ustrip.(bb), ustrip.(B))
+        end
+
+        # Scale invariance: renaming the units cannot change a dimensionless score.
+        # `A`'s cover is exact, so both scores sit at rounding level and only an
+        # absolute tolerance is meaningful here.
+        U2 = (u"mm", u"km/hr", u"kN")
+        A2 = [uconvert(unit(1 / (u * v)), A[i, j]) for (i, u) in enumerate(U2), (j, v) in enumerate(U2)]
+        for ϕ in ϕs
+            @test cover_objective(ϕ, symcover(A2), A2) ≈ cover_objective(ϕ, a, A) atol = 1e-12
+        end
+    end
+
 end