diff --git a/Project.toml b/Project.toml
index b18d325..2b415ff 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FixedEffects"
 uuid = "c8885935-8500-56a7-9867-7708b20db0eb"
-version = "3.5.0"
+version = "3.5.1"
 
 [deps]
 GroupedArrays = "6407cd72-fade-4a84-8a1e-56e431fc1533"
diff --git a/benchmarks/benchmark.jl b/benchmarks/benchmark.jl
index 275e4a1..cc4d583 100644
--- a/benchmarks/benchmark.jl
+++ b/benchmarks/benchmark.jl
@@ -2,8 +2,8 @@
 # Benchmarks for solve_residuals! across backends (:cpu, :CUDA, :Metal).
 # Problems are defined once; every available backend runs the same
 # specifications. Each specification residualizes 6 columns through one
-# solver, which is how FixedEffectModels calls this package (first run
-# includes compilation).
+# solver at tol = 1e-6, which is how FixedEffectModels calls this package by
+# default (first run includes compilation).
 # Usage: julia --project -t auto benchmarks/benchmark.jl
 ##############################################################################
 
@@ -64,33 +64,33 @@ println("="^60)
 
 feM = AbstractFixedEffectSolver{Float64}(fes_simple, uweights(length(cols_simple[1])), Val{:cpu})
 println("Simple (N=10M, 100k×100), 6 columns, first run:")
-@time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false)   # ~1.1 s
+@time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false, tol = 1e-6)   # ~0.73 s
 println("Simple (N=10M, 100k×100), 6 columns, second run:")
-@time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false)   # ~0.8 s
+@time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false, tol = 1e-6)   # ~0.60 s
 
 feM = AbstractFixedEffectSolver{Float64}(fes_hard, uweights(N), Val{:cpu})
 println("Hard (N=800k, 40k×5k), 6 columns, first run:")
-@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~13.6 s
+@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~3.0 s
 println("Hard (N=800k, 40k×5k), 6 columns, second run:")
-@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~13.8 s
+@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~3.1 s
 
 feM = AbstractFixedEffectSolver{Float64}(fes_hard_interact, uweights(N), Val{:cpu})
 println("Hard (N=800k, interacted), 6 columns, first run:")
-@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~22 s
+@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~6.1 s
 println("Hard (N=800k, interacted), 6 columns, second run:")
-@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~26 s
+@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~5.4 s
 
 feM = AbstractFixedEffectSolver{Float64}(fes_hard_3way, uweights(N), Val{:cpu})
 println("Hard 3-way (N=800k, 40k×5k×36), 6 columns, first run:")
-@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~20.6 s
+@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~4.1 s
 println("Hard 3-way (N=800k, 40k×5k×36), 6 columns, second run:")
-@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~19.9 s
+@time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~4.1 s
 
 feM = AbstractFixedEffectSolver{Float64}(fes_large, uweights(N_large), Val{:cpu})
 println("Large (N=10M, 500k×50k, sorted layout), 6 columns, first run:")
-@time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false)   # ~4.4 s
+@time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false, tol = 1e-6)   # ~2.7 s
 println("Large (N=10M, 500k×50k, sorted layout), 6 columns, second run:")
-@time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false)   # ~4.1 s
+@time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false, tol = 1e-6)   # ~2.7 s
 
 ##############################################################################
 # Metal
@@ -103,33 +103,33 @@ if isdefined(Main, :Metal) && Metal.functional()
 
     feM = AbstractFixedEffectSolver{Float32}(fes_simple, uweights(length(cols_simple[1])), Val{:Metal})
     println("Simple (N=10M, 100k×100), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false)   # ~9.9 s
+    @time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false, tol = 1e-6)   # ~4.3 s
     println("Simple (N=10M, 100k×100), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false)   # ~0.38 s
+    @time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false, tol = 1e-6)   # ~0.34 s
 
     feM = AbstractFixedEffectSolver{Float32}(fes_hard, uweights(N), Val{:Metal})
     println("Hard (N=800k, 40k×5k), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~3.5 s
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~2.8 s
     println("Hard (N=800k, 40k×5k), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~0.75 s
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~2.5 s
 
     feM = AbstractFixedEffectSolver{Float32}(fes_hard_interact, uweights(N), Val{:Metal})
     println("Hard (N=800k, interacted), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~3.2 s
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~4.2 s
     println("Hard (N=800k, interacted), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~1.3 s
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~3.9 s
 
     feM = AbstractFixedEffectSolver{Float32}(fes_hard_3way, uweights(N), Val{:Metal})
     println("Hard 3-way (N=800k, 40k×5k×36), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~1.6 s
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~3.6 s
     println("Hard 3-way (N=800k, 40k×5k×36), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)   # ~0.95 s
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)   # ~3.5 s
 
     feM = AbstractFixedEffectSolver{Float32}(fes_large, uweights(N_large), Val{:Metal})
     println("Large (N=10M, 500k×50k, sorted layout), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false)   # ~2.8 s
+    @time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false, tol = 1e-6)   # ~1.5 s
     println("Large (N=10M, 500k×50k, sorted layout), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false)   # ~1.0 s
+    @time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false, tol = 1e-6)   # ~1.3 s
 end
 
 ##############################################################################
@@ -143,31 +143,31 @@ if isdefined(Main, :CUDA) && CUDA.functional()
 
     feM = AbstractFixedEffectSolver{Float32}(fes_simple, uweights(length(cols_simple[1])), Val{:CUDA})
     println("Simple (N=10M, 100k×100), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false, tol = 1e-6)
     println("Simple (N=10M, 100k×100), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_simple], feM; progress_bar = false, tol = 1e-6)
 
     feM = AbstractFixedEffectSolver{Float32}(fes_hard, uweights(N), Val{:CUDA})
     println("Hard (N=800k, 40k×5k), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)
     println("Hard (N=800k, 40k×5k), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)
 
     feM = AbstractFixedEffectSolver{Float32}(fes_hard_interact, uweights(N), Val{:CUDA})
     println("Hard (N=800k, interacted), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)
     println("Hard (N=800k, interacted), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)
 
     feM = AbstractFixedEffectSolver{Float32}(fes_hard_3way, uweights(N), Val{:CUDA})
     println("Hard 3-way (N=800k, 40k×5k×36), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)
     println("Hard 3-way (N=800k, 40k×5k×36), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_hard], feM; progress_bar = false, tol = 1e-6)
 
     feM = AbstractFixedEffectSolver{Float32}(fes_large, uweights(N_large), Val{:CUDA})
     println("Large (N=10M, 500k×50k, sorted layout), 6 columns, first run:")
-    @time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false, tol = 1e-6)
     println("Large (N=10M, 500k×50k, sorted layout), 6 columns, second run:")
-    @time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false)
+    @time solve_residuals!([copy(c) for c in cols_large], feM; progress_bar = false, tol = 1e-6)
 end
diff --git a/ext/CUDAExt.jl b/ext/CUDAExt.jl
index ded5a67..f8b314d 100644
--- a/ext/CUDAExt.jl
+++ b/ext/CUDAExt.jl
@@ -1,7 +1,6 @@
 module CUDAExt
 using FixedEffects, CUDA
 using FixedEffects: FixedEffectCoefficients, AbstractWeights, UnitWeights, LinearAlgebra, Adjoint, mul!, rmul!, AbstractFixedEffectLinearMap, copy_internal!, AbsorptionPlan, AbsorbedBlock, _group_permutation, block_width
-CUDA.allowscalar(false)
 
 ##############################################################################
 ##
diff --git a/ext/MetalExt.jl b/ext/MetalExt.jl
index 25119a6..7541183 100644
--- a/ext/MetalExt.jl
+++ b/ext/MetalExt.jl
@@ -1,7 +1,6 @@
 module MetalExt
 using FixedEffects, Metal
 using FixedEffects: FixedEffectCoefficients, AbstractWeights, UnitWeights, LinearAlgebra, Adjoint, mul!, rmul!, AbstractFixedEffectLinearMap, copy_internal!, AbsorptionPlan, AbsorbedBlock, _group_permutation, block_width
-Metal.allowscalar(false)
 
 ##############################################################################
 ##
@@ -29,11 +28,6 @@ Metal.allowscalar(false)
 _mtl(T::Type, w::UnitWeights) = Metal.ones(T, length(w))
 _mtl(T::Type, w::AbstractVector) = MtlVector{T}(convert(Vector{T}, w))
 
-function _metal_threadgroup_width()
-	width = Int(device().maxThreadsPerThreadgroup.width)
-	return prevpow(2, width)
-end
-
 # Largest power-of-two threadgroup width the compiled pipeline allows: kernels
 # with heavy register use can cap below the device width.
 _pipeline_width(kernel) = prevpow(2, Int(kernel.pipeline.maxTotalThreadsPerThreadgroup))
@@ -52,12 +46,13 @@ function _sized_kernel(build, nthreads::Int)
 	return kernel, nthreads
 end
 
-# Per-block plan for the adjoint gather (A'u), chosen once at construction:
-# bucketize (one threadgroup per group) for low cardinality, else atomic adds.
+# Per-block plan for the adjoint gather (A'u), chosen once at construction by
+# _gather_strategy: atomic adds in row order, or one threadgroup per group.
 struct AtomicGather end
 struct BucketGather{V<:AbstractVector}
 	perm::V
 	offsets::V
+	width::Int
 end
 
 mutable struct FixedEffectLinearMapMetal{T,P<:AbsorptionPlan} <: AbstractFixedEffectLinearMap{T}
@@ -77,13 +72,30 @@ function FixedEffectLinearMapMetal{T}(fes::Vector{<:FixedEffect}, weights::Abstr
 	return FixedEffectLinearMapMetal{T,typeof(plan)}(fes, plan, gathers)
 end
 
+# Atomic adds are the fastest gather while threads rarely update the same
+# coefficient at once, but they become orders of magnitude slower when they do,
+# which happens with a large group or with rows clustered by group. They are
+# used only when neither applies. Otherwise each group gets a threadgroup, whose
+# width trades off idle threads (nlevels * width) against the time spent in the
+# largest group (nmax / width). Widths are powers of 4 to limit the number of
+# kernels to compile. The constants were calibrated on an M1 Max.
 function _gather_strategy(refs::AbstractVector{<:Integer}, nlevels::Int)
-	if nlevels < min(100_000, div(length(refs), 16))
-		_, offsets, perm = _group_permutation(refs, nlevels)
-		return BucketGather(MtlVector{Int}(perm), MtlVector{Int}(offsets))
-	else
+	counts = zeros(Int, nlevels)
+	nadjacent = 0
+	previous = zero(eltype(refs))
+	@inbounds for g in refs
+		counts[g] += 1
+		nadjacent += g == previous
+		previous = g
+	end
+	nmax = maximum(counts; init = 0)
+	if nmax <= div(length(refs), 20_000) && 2 * nadjacent <= length(refs)
 		return AtomicGather()
 	end
+	device_width = prevpow(2, Int(device().maxThreadsPerThreadgroup.width))
+	width = prevpow(4, clamp(round(Int, 90 * sqrt(nmax / nlevels)), 16, device_width))
+	_, offsets, perm = _group_permutation(refs, nlevels)
+	return BucketGather(MtlVector{Int}(perm), MtlVector{Int}(offsets), width)
 end
 
 function _mtl_plan(::Type{T}, fes::Vector{<:FixedEffect}, weights::AbstractWeights) where {T}
@@ -95,45 +107,85 @@ function _mtl_plan(::Type{T}, fes::Vector{<:FixedEffect}, weights::AbstractWeigh
 end
 
 ## 1b) FixedEffectLinearMapMetal mul!
+##
+## One kernel computes u ← c u + α Σ_blocks A_b v_b for all blocks at once and
+## accumulates ‖u‖² on the fly: one partial per threadgroup, summed pairwise on
+## the device — Metal has no Float64 atomics, and the pairwise sum avoids the
+## drift of chaining Float32 atomic adds (see the _norm2 note in
+## src/utils/lsmr.jl). Blocks are passed as tuples of device arrays so the
+## kernel unrolls across them. mul! discards the norm; bidiag_forward! uses it.
 
 ## Implement right multiplication
 function LinearAlgebra.mul!(y::MtlVector, fem::FixedEffectLinearMapMetal{T},
 		fecoefs::FixedEffectCoefficients, α::Number, β::Number) where {T}
-	if iszero(β)
-		fill!(y, zero(T))
-		β = one(β)
-	end
-	for (coef_block, block, qrows) in zip(fecoefs.x, fem.plan.blocks, fem.plan.qrows)
-		_scatter_block!(y, block.refs, qrows, coef_block, α, β)
-		β = one(β)
-	end
+	# β = 0 must ignore the contents of y, which may be NaN
+	iszero(β) && fill!(y, zero(T))
+	_scatter!(y, fem, fecoefs, α, β)
 	return y
 end
 
-function _scatter_block!(y::MtlVector, refs::MtlVector, qrows::MtlMatrix,
-		coef_block::MtlMatrix, α::Number, β::Number)
-	kernel = @metal launch=false scatter_block_kernel!(y, refs, qrows, coef_block, α, β, size(coef_block, 1))
-	nthreads = _pipeline_width(kernel)
-	kernel(y, refs, qrows, coef_block, α, β, size(coef_block, 1);
-		threads = nthreads, groups = cld(length(y), nthreads))
-	return y
+# u ← c u + α A v; returns the per-threadgroup partial sums of u².
+function _scatter!(u::MtlVector{T}, fem::FixedEffectLinearMapMetal{T},
+		v::FixedEffectCoefficients, α::Number, c::Number) where {T}
+	refss = map(block -> block.refs, Tuple(fem.plan.blocks))
+	qrowss = Tuple(fem.plan.qrows)
+	vs = Tuple(v.x)
+	# fixed 256-thread groups, safely below any pipeline cap for this kernel
+	nthreads = 256
+	nblocks = cld(length(u), nthreads)
+	partials = MtlVector{T}(undef, nblocks)
+	@metal threads=nthreads groups=nblocks scatter_kernel!(u, refss, qrowss, vs, T(α), T(c),
+		partials, Val(nthreads))
+	return partials
 end
 
-function scatter_block_kernel!(y, refs, qrows, coef_block, α, β, k)
+function scatter_kernel!(u, refss, qrowss, vs, α, c, partials, ::Val{NT}) where {NT}
+	T = eltype(u)
+	gid = Int(threadgroup_position_in_grid().x)
+	tid = Int(thread_position_in_threadgroup().x)
+	nt = Int(threads_per_threadgroup().x)
+	shared = Metal.MtlThreadGroupArray(T, NT)
+	acc = zero(T)
 	i = thread_position_in_grid_1d()
-	if i <= length(y)
+	if i <= length(u)
 		@inbounds begin
-			g = refs[i]
-			fit = zero(eltype(y))
-			for c in 1:k
-				fit += coef_block[c, g] * qrows[c, i]
-			end
-			y[i] = β * y[i] + α * fit
+			ui = c * u[i] + α * _device_fit(refss, qrowss, vs, i)
+			u[i] = ui
+			acc += ui * ui
+		end
+	end
+	@inbounds shared[tid] = acc
+	Metal.threadgroup_barrier(Metal.MemoryFlagThreadGroup)
+
+	offset = nt ÷ 2
+	while offset > 0
+		if tid <= offset
+			@inbounds shared[tid] += shared[tid + offset]
 		end
+		Metal.threadgroup_barrier(Metal.MemoryFlagThreadGroup)
+		offset ÷= 2
+	end
+
+	if tid == 1
+		@inbounds partials[gid] = shared[1]
 	end
 	return nothing
 end
 
+# Recursion over the block tuples, as in the CPU kernels.
+@inline _device_fit(::Tuple{}, ::Tuple{}, ::Tuple{}, i) = false
+@inline function _device_fit(refss::Tuple, qrowss::Tuple, vs::Tuple, i)
+	refs = first(refss)
+	qrows = first(qrowss)
+	vcoef = first(vs)
+	@inbounds gr = refs[i]
+	fit = zero(eltype(qrows))
+	for col in 1:size(qrows, 1)
+		@inbounds fit += vcoef[col, gr] * qrows[col, i]
+	end
+	return fit + _device_fit(Base.tail(refss), Base.tail(qrowss), Base.tail(vs), i)
+end
+
 ## 1c) FixedEffectLinearMapMetal mul!, Adjoint
 
 ## Implement left multiplication
@@ -151,7 +203,7 @@ end
 function _gather_block!(coef_block::MtlMatrix, refs::MtlVector, qrows::MtlMatrix,
 		y::MtlVector, α::Number, gather::BucketGather)
 	k = size(coef_block, 1)
-	kernel, nthreads = _sized_kernel(_metal_threadgroup_width()) do nt
+	kernel, nthreads = _sized_kernel(gather.width) do nt
 		@metal launch=false gather_block_kernel_bin!(coef_block, α, y, qrows, gather.perm, gather.offsets, Val(nt), k)
 	end
 	kernel(coef_block, α, y, qrows, gather.perm, gather.offsets, Val(nthreads), k;
@@ -226,27 +278,15 @@ end
 ##
 ## 1d) Fused bidiagonalization step (single-pass LSMR iterations)
 ##
-## One kernel computes u ← Σ_blocks A_b v_b + c u and accumulates ‖u‖² on the
-## fly: one partial per threadgroup, summed pairwise on the device — Metal has
-## no Float64 atomics, and the pairwise sum avoids the drift of chaining
-## Float32 atomic adds (see the _norm2 note in src/utils/lsmr.jl). The gathers
-## then run through their per-block strategies on the raw u. Blocks are passed
-## as tuples of device arrays so the kernel unrolls across them.
+## The scatter kernel of 1b computes u ← A v + c u and ‖u‖² in one pass over
+## the rows; the gathers then run through their per-block strategies on the
+## raw u.
 ##
 ##############################################################################
 
 function FixedEffects.bidiag_forward!(u::MtlVector{T}, g::FixedEffectCoefficients,
 		fem::FixedEffectLinearMapMetal{T}, v::FixedEffectCoefficients, c::Number) where {T}
-	blocks = Tuple(fem.plan.blocks)
-	refss = map(block -> block.refs, blocks)
-	qrowss = Tuple(fem.plan.qrows)
-	vs = Tuple(v.x)
-	# fixed 256-thread groups, safely below any pipeline cap for this kernel
-	nthreads = 256
-	nblocks = cld(length(u), nthreads)
-	partials = MtlVector{T}(undef, nblocks)
-	@metal threads=nthreads groups=nblocks bidiag_scatter_kernel!(u, refss, qrowss, vs, T(c),
-		partials, Val(nthreads))
+	partials = _scatter!(u, fem, v, one(T), c)
 	fill!(g, zero(T))
 	for (coef_block, block, qrows, gather) in zip(g.x, fem.plan.blocks, fem.plan.qrows, fem.gathers)
 		_gather_block!(coef_block, block.refs, qrows, u, one(T), gather)
@@ -256,53 +296,6 @@ function FixedEffects.bidiag_forward!(u::MtlVector{T}, g::FixedEffectCoefficient
 	return sqrt(sum(partials))
 end
 
-function bidiag_scatter_kernel!(u, refss, qrowss, vs, c, partials, ::Val{NT}) where {NT}
-	T = eltype(u)
-	gid = Int(threadgroup_position_in_grid().x)
-	tid = Int(thread_position_in_threadgroup().x)
-	nt = Int(threads_per_threadgroup().x)
-	shared = Metal.MtlThreadGroupArray(T, NT)
-	acc = zero(T)
-	i = thread_position_in_grid_1d()
-	if i <= length(u)
-		@inbounds begin
-			ui = c * u[i] + _device_fit(refss, qrowss, vs, i)
-			u[i] = ui
-			acc += ui * ui
-		end
-	end
-	@inbounds shared[tid] = acc
-	Metal.threadgroup_barrier(Metal.MemoryFlagThreadGroup)
-
-	offset = nt ÷ 2
-	while offset > 0
-		if tid <= offset
-			@inbounds shared[tid] += shared[tid + offset]
-		end
-		Metal.threadgroup_barrier(Metal.MemoryFlagThreadGroup)
-		offset ÷= 2
-	end
-
-	if tid == 1
-		@inbounds partials[gid] = shared[1]
-	end
-	return nothing
-end
-
-# Recursion over the block tuples, as in the CPU kernels.
-@inline _device_fit(::Tuple{}, ::Tuple{}, ::Tuple{}, i) = false
-@inline function _device_fit(refss::Tuple, qrowss::Tuple, vs::Tuple, i)
-	refs = first(refss)
-	qrows = first(qrowss)
-	vcoef = first(vs)
-	@inbounds gr = refs[i]
-	fit = zero(eltype(qrows))
-	for col in 1:size(qrows, 1)
-		@inbounds fit += vcoef[col, gr] * qrows[col, i]
-	end
-	return fit + _device_fit(Base.tail(refss), Base.tail(qrowss), Base.tail(vs), i)
-end
-
 ##############################################################################
 ##
 ## 2. FixedEffectSolverMetal
diff --git a/src/AbstractFixedEffectSolver.jl b/src/AbstractFixedEffectSolver.jl
index e94909f..d4de7a0 100644
--- a/src/AbstractFixedEffectSolver.jl
+++ b/src/AbstractFixedEffectSolver.jl
@@ -50,7 +50,7 @@ function solve_residuals!(y::AbstractVector{<: Real}, fes::AbstractVector{<: Fix
 	solve_residuals!(y, feM; maxiter = maxiter, tol = tol)
 end
 
-function solve_residuals!(r::AbstractVector{<:Real}, feM::AbstractFixedEffectSolver{T}; tol::Real = sqrt(eps(T)), maxiter::Integer = 100_000) where {T}
+function solve_residuals!(r::AbstractVector{<:Real}, feM::AbstractFixedEffectSolver{T}; tol::Real = T === Float64 ? 1e-8 : 1e-6, maxiter::Integer = 100_000) where {T}
 	maxiter >= 0 || throw(ArgumentError("maxiter must be non-negative"))
 	# One cannot copy view of Vector (r) on GPU, so first collect the vector
 	copy_internal!(feM, :r, r)
@@ -155,7 +155,7 @@ function solve_coefficients!(y::AbstractVector{<: Number}, fes::AbstractVector{<
 	solve_coefficients!(y, feM; maxiter = maxiter, tol = tol)
 end
 
-function solve_coefficients!(r::AbstractVector, feM::AbstractFixedEffectSolver{T}; tol::Real = sqrt(eps(T)), maxiter::Integer = 100_000) where {T}
+function solve_coefficients!(r::AbstractVector, feM::AbstractFixedEffectSolver{T}; tol::Real = T === Float64 ? 1e-8 : 1e-6, maxiter::Integer = 100_000) where {T}
 	maxiter >= 0 || throw(ArgumentError("maxiter must be non-negative"))
 	# One cannot copy view of Vector (r) on GPU, so first collect the vector
 	copy_internal!(feM, :b, r)
diff --git a/test/solve.jl b/test/solve.jl
index 8738bc9..c8cd780 100644
--- a/test/solve.jl
+++ b/test/solve.jl
@@ -43,6 +43,19 @@ end
 	@test_throws ArgumentError solve_coefficients!(deepcopy(x), fes; maxiter = -1)
 end
 
+@testset "default tol is the same for both entry points" begin
+	n = 2_000
+	pid = [mod(i^2, 97) + 1 for i in 1:n]
+	fid = [div(pid[i], 8) + 1 + mod(i, 3) for i in 1:n]
+	fes_tol = [FixedEffect(pid), FixedEffect(fid)]
+	y_tol = [sin(i) for i in 1:n]
+	for (T, double_precision) in ((Float64, true), (Float32, false))
+		feM_tol = FixedEffects.AbstractFixedEffectSolver{T}(fes_tol, uweights(T, n), Val{:cpu})
+		@test solve_residuals!(copy(y_tol), feM_tol)[2] == solve_residuals!(copy(y_tol), fes_tol; double_precision = double_precision)[2]
+		@test solve_coefficients!(copy(y_tol), feM_tol)[2] == solve_coefficients!(copy(y_tol), fes_tol; double_precision = double_precision)[2]
+	end
+end
+
 # PooledArrays
 (r, iter, conv) = solve_residuals!(deepcopy(x), [FixedEffect(PooledArray(p1)), FixedEffect(PooledArray(p2))])
 @test r ≈ r_ols
@@ -207,6 +220,29 @@ end
 	end
 end
 
+if Metal.functional()
+	@testset "Metal gather strategy" begin
+		ext = Base.get_extension(FixedEffects, :MetalExt)
+		n = 100_000
+		scattered = mod1.(7919 .* (1:n), 20_000)
+		clustered = sort(scattered)
+		skewed = copy(scattered)
+		skewed[1:2:n] .= 1
+		few = mod1.(1:n, 4)
+		@test ext._gather_strategy(scattered, 20_000) isa ext.AtomicGather
+		@test ext._gather_strategy(clustered, 20_000).width == 16
+		@test ext._gather_strategy(skewed, 20_000).width == 64
+		@test ext._gather_strategy(few, 4).width > 64
+		y_gather = sin.(1:n)
+		for refs in (scattered, clustered, skewed, few)
+			fes_gather = [FixedEffect(refs), FixedEffect(mod1.(1:n, 7))]
+			r_cpu = solve_residuals!(copy(y_gather), fes_gather; double_precision = false)[1]
+			r_metal = solve_residuals!(copy(y_gather), fes_gather; method = :Metal, double_precision = false)[1]
+			@test r_metal ≈ r_cpu rtol = 1e-3
+		end
+	end
+end
+
 
 fe = FixedEffect([1, 2])
 @test_throws "FixedEffects must have the same length as y" ỹ = solve_residuals!(ones(100), [fe])