diff --git a/bench/fact_fused.jl b/bench/fact_fused.jl
index 58a1cc6..97c3c5d 100644
--- a/bench/fact_fused.jl
+++ b/bench/fact_fused.jl
@@ -411,6 +411,7 @@ function _sorted_subtree_launch(N::SDS.Numeric{T}, S, nzval, list, ::Val{LB}, ::
     k!(N.factor, N.stack, N.info, nzval, S.amap, S.amap_ptr, S.amap_src, list, bm, S.subtree_ptr,
        S.subtree_nodes, S.front_ptr, S.front_nrows, S.front_ncols, S.cb_ptr, S.local_front, S.local_cb,
        S.child_ptr, S.child_list, S.relind_ptr, S.relind,
+       real(T)(SDS.resolved_pivot_epsilon(SDS.Options(), real(T))),   # production default (#127)
        Val((LB - SDS.SUBTREE_LOCAL_RESERVE) ÷ sizeof(T)), Val(WG); ndrange = WG * length(list) * bm.nact)
     return nothing
 end
diff --git a/bench/fact_split.jl b/bench/fact_split.jl
index c4b7d38..cc229cb 100644
--- a/bench/fact_split.jl
+++ b/bench/fact_split.jl
@@ -248,6 +248,7 @@ function _sorted_subtree_launch2(N::SDS.Numeric{T}, S, nzval, list, bm, k!, ::Va
     k!(N.factor, N.stack, N.info, nzval, S.amap, S.amap_ptr, S.amap_src, list, bm, S.subtree_ptr,
        S.subtree_nodes, S.front_ptr, S.front_nrows, S.front_ncols, S.cb_ptr, S.local_front, S.local_cb,
        S.child_ptr, S.child_list, S.relind_ptr, S.relind,
+       real(T)(SDS.resolved_pivot_epsilon(SDS.Options(), real(T))),   # production default (#127)
        Val((LB - SDS.SUBTREE_LOCAL_RESERVE) ÷ sizeof(T)), Val(WG); ndrange = WG * length(list) * bm.nact)
     return nothing
 end
diff --git a/src/dense/fallback/potrf.jl b/src/dense/fallback/potrf.jl
index 71a2bb1..cfd3a62 100644
--- a/src/dense/fallback/potrf.jl
+++ b/src/dense/fallback/potrf.jl
@@ -113,7 +113,7 @@ function ka_chol_diag_info!(info::AbstractVector{Int32}, L::AbstractArray)
     return info
 end
 
-@kernel function _ka_chol_check_kernel!(info, A, n, idx)
+@kernel function _ka_chol_check_kernel!(info, A, n, idx, dmin)
     li = @index(Local, Linear)
     b = @index(Group, Linear)
     if li == 1
@@ -121,7 +121,11 @@ end
             st = Int32(0)
             for j in 1:n
                 d = _get(A, j, j, b)
-                if !(isfinite(real(d)) && real(d) > 0 && iszero(imag(d)))
+                # dmin > 0: the Cholesky pivot-epsilon floor (issue #127) — potrf accepts
+                # any positive pivot, so a pivot that became tiny DURING elimination is
+                # only visible here, as a factor diagonal at or below sqrt(epsilon)
+                ok = dmin > zero(dmin) ? real(d) >= dmin : real(d) > zero(real(d))
+                if !(isfinite(real(d)) && ok && iszero(imag(d)))
                     st = Int32(j)
                     break
                 end
@@ -140,10 +144,10 @@ the first `j` whose diagonal entry `L[j, j]` is not a finite positive real
 (still 0 if there is none); a nonzero `info[idx]` is kept. Asynchronous (see
 [`potrf_info!`](@ref)).
 """
-function ka_chol_check_info!(info::AbstractVector{Int32}, idx::Integer, L::AbstractMatrix)
+function ka_chol_check_info!(info::AbstractVector{Int32}, idx::Integer, L::AbstractMatrix, dmin::Real = 0.0)
     1 <= idx <= length(info) || throw(DimensionMismatch("info index $idx outside 1:$(length(info))"))
     kernel! = _ka_chol_check_kernel!(KernelAbstractions.get_backend(L), 1)
-    kernel!(info, L, size(L, 1), Int(idx); ndrange = 1)
+    kernel!(info, L, size(L, 1), Int(idx), real(eltype(L))(dmin); ndrange = 1)
     return info
 end
 
@@ -153,12 +157,13 @@ end
 Batched form: the check of `ka_chol_check_info!(info, offset + b, L[:, :, b])`
 for every member `b` of the 3-D strided batch `L`, one launch.
 """
-function ka_chol_check_info!(info::AbstractVector{Int32}, offset::Integer, L::AbstractArray{<:Any, 3})
+function ka_chol_check_info!(info::AbstractVector{Int32}, offset::Integer, L::AbstractArray{<:Any, 3},
+                             dmin::Real = 0.0)
     nb = size(L, 3)
     0 <= offset && offset + nb <= length(info) ||
         throw(DimensionMismatch("info has length $(length(info)) < offset $offset + batch count $nb"))
     nb == 0 && return info
     kernel! = _ka_chol_check_kernel!(KernelAbstractions.get_backend(L), 1)
-    kernel!(info, L, size(L, 1), Int(offset) + 1; ndrange = nb)
+    kernel!(info, L, size(L, 1), Int(offset) + 1, real(eltype(L))(dmin); ndrange = nb)
     return info
 end
diff --git a/src/dense/interface.jl b/src/dense/interface.jl
index c53543d..21ad7c3 100644
--- a/src/dense/interface.jl
+++ b/src/dense/interface.jl
@@ -287,16 +287,19 @@ function potrf_info!(uplo, A::AbstractMatrix, info::AbstractVector{Int32}, idx::
 end
 
 # `potrf_info!` with an already resolved implementation `p` (no capability lookup)
-function _potrf_info_impl!(p::Symbol, ul::Char, A::AbstractMatrix, info::AbstractVector{Int32}, idx::Integer)
+function _potrf_info_impl!(p::Symbol, ul::Char, A::AbstractMatrix, info::AbstractVector{Int32}, idx::Integer,
+                           dmin::Real = 0.0)
     if p === :vendor
         vendor_potrf_info!(ul, A, info, idx)
-        ka_chol_check_info!(info, idx, A)
+        ka_chol_check_info!(info, idx, A, dmin)
     elseif p === :generic
         cholesky!(Hermitian(A, ul == 'L' ? :L : :U), NoPivot(); check = false)
         fill!(view(info, idx:idx), Int32(0))
-        ka_chol_check_info!(info, idx, A)
+        ka_chol_check_info!(info, idx, A, dmin)
     else
         ka_potrf!(ul, A, info; offset = idx - 1)
+        # ka_potrf validates its own pivots at > 0; the epsilon floor needs the pass too
+        dmin > 0 && ka_chol_check_info!(info, idx, A, dmin)
     end
     return info
 end
diff --git a/src/fused.jl b/src/fused.jl
index f692394..a0c6a69 100644
--- a/src/fused.jl
+++ b/src/fused.jl
@@ -992,7 +992,7 @@ end
                                     cb_done, @Const(snparent), @Const(front_ptr), @Const(front_nrows),
                                     @Const(front_ncols), @Const(cb_ptr), @Const(child_ptr),
                                     @Const(child_list), @Const(relind_ptr), @Const(relind), base,
-                                    @Const(wbase), wupd, wpdone2, wtrsml, wcbl, scap, eac, d, @Const(super_ptr),
+                                    @Const(wbase), wupd, wpdone2, wtrsml, wcbl, scap, eac, epsc, d, @Const(super_ptr),
                                     ::Val{LDL}, ::Val{WG}) where {LDL, WG}
     @uniform TT = eltype(factor)
     @uniform RT = real(eltype(factor))
@@ -1103,13 +1103,13 @@ end
                 if LDL
                     _front_ldl_update!(sh, st, piv, j, v, li, front_ncols, d, g0, Val(WG))
                 else
-                    SDS._front_chol_update!(sh, st, piv, j, v, li, front_ncols, Val(64), Val(WG))
+                    SDS._front_chol_update!(sh, st, piv, j, v, li, front_ncols, epsc, Val(64), Val(WG))
                 end
                 @synchronize
                 if LDL
                     _front_ldl_scale!(sh, st, piv, j, v, li, front_ncols, Val(WG))
                 else
-                    SDS._front_chol_scale!(sh, st, piv, j, v, li, front_ncols, Val(64), Val(WG))
+                    SDS._front_chol_scale!(sh, st, piv, j, v, li, front_ncols, epsc, Val(64), Val(WG))
                 end
                 @synchronize
             end
@@ -1365,8 +1365,10 @@ end
                             piv[1] = one(RT)
                         end
                     else
-                        if st[1] == Int32(0) && dj > zero(RT)
-                            piv[1] = sqrt(dj)
+                        # same pivot-epsilon semantics as the narrow kernels (#127):
+                        # [0, epsc] perturbs to epsc, negative/NaN fails the front
+                        if st[1] == Int32(0) && (dj > epsc || (epsc > zero(RT) && dj >= zero(RT)))
+                            piv[1] = sqrt(max(dj, epsc))
                             sh[SDS._packed(j, j, 64)] = piv[1]
                         elseif st[1] == Int32(0)
                             st[1] = Int32(c0 + j)
@@ -1765,13 +1767,13 @@ end
 end
 
 # subtree launcher with a custom (sorted) list, local memory sized by eltype
-function launch_subtrees_sorted!(N, S, nzval, list, lb, backend, ldl::Bool = false)
+function launch_subtrees_sorted!(N, S, nzval, list, lb, backend, ldl::Bool = false, epsc::Real = 0.0)
     SDS._with_local_bytes(Int(lb)) do LBv
-        _sorted_subtree_launch(N, S, nzval, list, Val(ldl), LBv, Val(SDS.SUBTREE_WORKGROUP), backend)
+        _sorted_subtree_launch(N, S, nzval, list, Val(ldl), LBv, Val(SDS.SUBTREE_WORKGROUP), backend, epsc)
     end
 end
 function _sorted_subtree_launch(N::SDS.Numeric{T}, S, nzval, list, ::Val{LDL}, ::Val{LB}, ::Val{WG},
-                                backend) where {T, LDL, LB, WG}
+                                backend, epsc::Real = 0.0) where {T, LDL, LB, WG}
     bm = SDS.batch_map(N; first = 1)
     if LDL
         # the schedule sized the "S" budget classes with the LDLT reserve: match it
@@ -1785,7 +1787,7 @@ function _sorted_subtree_launch(N::SDS.Numeric{T}, S, nzval, list, ::Val{LDL}, :
         k! = SDS.subtree_cholesky_kernel!(backend, WG)
         k!(N.factor, N.stack, N.info, nzval, S.amap, S.amap_ptr, S.amap_src, list, bm, S.subtree_ptr,
            S.subtree_nodes, S.front_ptr, S.front_nrows, S.front_ncols, S.cb_ptr, S.local_front, S.local_cb,
-           S.child_ptr, S.child_list, S.relind_ptr, S.relind,
+           S.child_ptr, S.child_list, S.relind_ptr, S.relind, real(T)(epsc),
            Val((LB - SDS.SUBTREE_LOCAL_RESERVE) ÷ sizeof(T)), Val(WG); ndrange = WG * length(list) * bm.nact)
     end
     return nothing
@@ -2021,6 +2023,9 @@ function build_fused_plan(s; ldl::Bool = false, inkernel_wides::Union{Nothing, B
             wcb_b = _dvec(backend, wcb_b), wupd = KA.zeros(backend, Int32, max(nwslots, 1)),
             wtl = KA.zeros(backend, Int32, max(nwslots, 1)), wpd = KA.zeros(backend, Int32, max(nwslots, 1)),
             wcb = KA.zeros(backend, Int32, ns), wgf, scap = Int32(scap), eac = Int32(ea_chunk), ldl,
+            # captured at plan build: a later setparam!("pivot_epsilon") changes the stock
+            # path but not an already-built fused plan (rebuild the plan to change it)
+            chol_eps = ldl ? 0.0 : Float64(SDS.resolved_pivot_epsilon(s.options, real(T))),
             ldl_potrf = ldl && wides_potrf,
             sptr_h = Array(S.super_ptr),
             ldlw = KA.zeros(backend, T,
@@ -2050,7 +2055,7 @@ function refact_fused!(s, nzval, fp)
     fill!(fp.wpd, Int32(0))
     copyto!(fp.wcb, fp.wcb_b)
     for k in eachindex(plan.sub_first)
-        launch_subtrees_sorted!(N, S, nzval, fp.ssub[k], plan.sub_local[k], backend, fp.ldl)
+        launch_subtrees_sorted!(N, S, nzval, fp.ssub[k], plan.sub_local[k], backend, fp.ldl, fp.chol_eps)
     end
     nodes = S.schedule.group_nodes
     fill!(fp.pcb, zero(eltype(fp.pcb)))
@@ -2070,7 +2075,8 @@ function refact_fused!(s, nzval, fp)
                                              S.snparent, S.front_ptr, S.front_nrows, S.front_ncols,
                                              S.cb_ptr, S.child_ptr, S.child_list, S.relind_ptr, S.relind,
                                              Int32(base), fp.wbase, fp.wupd, fp.wpd, fp.wtl, fp.wcb,
-                                             fp.scap, fp.eac, N.d, S.super_ptr, Val(fp.ldl), Val(fp.wgf);
+                                             fp.scap, fp.eac, real(eltype(N.factor))(fp.chol_eps), N.d,
+                                             S.super_ptr, Val(fp.ldl), Val(fp.wgf);
                                              ndrange = fp.wgf * count)
         for v in gwides                              # host-wide fronts at the segment boundary
             if fp.ldl && fp.ldl_potrf
@@ -2079,7 +2085,7 @@ function refact_fused!(s, nzval, fp)
                 _factor_panel_c_ldl!(N, S, fp, v, p, backend)
             else
                 SDS._factor_panel_c!(N.factor, fp.pcb, N.work, N.info, v, S.layout.panel_ptr[v],
-                                     S.schedule.rows[v], S.schedule.width[v], fp.cb2_h[v], p)
+                                     S.schedule.rows[v], S.schedule.width[v], fp.cb2_h[v], p, fp.chol_eps)
             end
             signal_kernel!(backend, 32)(fp.cbd, fp.arrived, S.snparent, v; ndrange = 32)
         end
diff --git a/src/numeric/factorize.jl b/src/numeric/factorize.jl
index f784977..52b3e3c 100644
--- a/src/numeric/factorize.jl
+++ b/src/numeric/factorize.jl
@@ -109,11 +109,11 @@ end
 
 # regime-C dense kernels of front `s` (assembled panel + contribution block) of every active batch member,
 # `p` from `_front_impls`
-function _factor_front_c!(N::Numeric{T}, S::Symbolic, s::Int, p::NamedTuple) where {T}
+function _factor_front_c!(N::Numeric{T}, S::Symbolic, s::Int, p::NamedTuple, epsc::Real = 0.0) where {T}
     L, sc = S.layout, S.schedule
     if N.nbatch == 1
         _factor_panel_c!(N.factor, N.stack, N.work, N.info, s, L.panel_ptr[s], sc.rows[s], sc.width[s], L.cb_ptr[s],
-                         p)
+                         p, epsc)
         return nothing
     end
     plan = N.plan
@@ -125,10 +125,19 @@ function _factor_front_c!(N::Numeric{T}, S::Symbolic, s::Int, p::NamedTuple) whe
             for j in j0:j1
                 k = Int(plan.members_host[j])
                 _factor_panel_c!(N.factor, _mview(N.stack, k, nb), N.work, N.info, (s - 1) * nb + k,
-                                 panel_offset(L.panel_ptr, s, k, nb), sc.rows[s], sc.width[s], L.cb_ptr[s], p)
+                                 panel_offset(L.panel_ptr, s, k, nb), sc.rows[s], sc.width[s], L.cb_ptr[s], p, epsc)
             end
         else
-            _factor_panels_c_batched!(N, S, s, j0, j1 - j0 + 1, p)
+            if epsc > 0
+                backend = KernelAbstractions.get_backend(N.factor)
+                for j in j0:j1
+                    k = Int(plan.members_host[j])
+                    _chol_diag_clamp_kernel!(backend, 64)(N.factor, panel_offset(L.panel_ptr, s, k, nb),
+                                                          sc.rows[s], sc.width[s], real(T)(epsc);
+                                                          ndrange = sc.width[s])
+                end
+            end
+            _factor_panels_c_batched!(N, S, s, j0, j1 - j0 + 1, p, epsc)
         end
     end
     return nothing
@@ -138,7 +147,7 @@ end
 # are one `f×w×count` strided batch (interleaved layout), factored by strided-batched `potrf`, `trsm` and `gemm`
 # (`F21 F21ᴴ` into `count` blocks of the workspace), then packed into the members' contribution blocks
 function _factor_panels_c_batched!(N::Numeric{T}, S::Symbolic, s::Int, j0::Int, count::Int,
-                                   p::NamedTuple) where {T}
+                                   p::NamedTuple, epsc::Real = 0.0) where {T}
     L, sc = S.layout, S.schedule
     nb = N.nbatch
     f, w = sc.rows[s], sc.width[s]
@@ -149,12 +158,14 @@ function _factor_panels_c_batched!(N::Numeric{T}, S::Symbolic, s::Int, j0::Int,
     F11 = view(P3, 1:w, 1:w, :)
     ioff = (s - 1) * nb + k0 - 1                              # info of front s, member k0 (interleaved)
     ptrs = N.plan.vendor_ptrs[s]                              # (F11, F21) member pointers, or nothing
+    dmin = sqrt(real(T)(epsc))                                # epsilon floor (#127); 0 keeps today's check
     if p.bpotrf === :vendor && ptrs !== nothing
         vendor_potrf_batched_ptrs!('L', w, view(ptrs[1], k0:(k0 + count - 1)), f,
                                    view(N.info, (ioff + 1):(ioff + count)), count)
-        ka_chol_check_info!(N.info, ioff, F11)
+        ka_chol_check_info!(N.info, ioff, F11, dmin)
     else
         _potrf_batched_impl!(p.bpotrf, F11, N.info, ioff)
+        epsc > 0 && ka_chol_check_info!(N.info, ioff, F11, dmin)
     end
     m > 0 || return nothing
     F21 = view(P3, (w + 1):f, 1:w, :)
@@ -175,13 +186,36 @@ end
 # `potrf`/`trsm`/`syrk` (`herk`) on the `f×w` panel at `factor[p0]`; the `syrk` result goes to the
 # `m×m` workspace `work[1:m^2]` and is added to the packed contribution block at `stack[c0]`
 # (`c0 = 0`: none); status into `info[s]` (also used by bench/front_bins.jl)
+# pivots in [0, epsc] of the still-unfactored panel diagonal are raised to epsc before
+# potrf (the cuDSS pivot_epsilon static perturbation on the vendor/host Cholesky path);
+# nonpositive diagonals are left for potrf to reject, as in the kernel paths
+@kernel function _chol_diag_clamp_kernel!(factor, p0::Int, f::Int, w::Int, epsc)
+    j = @index(Global)
+    @inbounds if j <= w
+        idx = p0 + (j - 1) * f + j - 1
+        d = real(factor[idx])
+        # [0, epsc]: the kernel routes perturb an exact zero too (cuDSS semantics)
+        (d >= zero(d)) & (d <= epsc) && (factor[idx] = epsc)
+    end
+end
+
 function _factor_panel_c!(factor::AbstractVector{T}, stack::AbstractVector{T}, work::AbstractVector{T},
                           info::AbstractVector{Int32}, s::Int, p0::Int, f::Int, w::Int, c0::Int,
-                          p::NamedTuple) where {T}
+                          p::NamedTuple, epsc::Real = 0.0) where {T}
     m = f - w
     P = reshape(view(factor, p0:(p0 + f * w - 1)), f, w)
     F11 = view(P, 1:w, 1:w)
-    _potrf_info_impl!(p.potrf, 'L', F11, info, s)
+    if epsc > 0
+        # same-queue ordering covers every potrf impl (the CPU backend already ran
+        # LAPACK straight after the assembly kernels before this PR): no host sync
+        _chol_diag_clamp_kernel!(KernelAbstractions.get_backend(factor), 64)(factor, p0, f, w, real(T)(epsc);
+                                                                             ndrange = w)
+    end
+    # the clamp handles input-diagonal tininess; a pivot that becomes tiny DURING the
+    # elimination is caught after the fact by the sqrt(epsilon) factor-diagonal floor
+    # inside _potrf_info_impl!, reported through info (the potrf route DETECTS where the
+    # kernel routes perturb -- either way MadNLP gets its regularization signal, #127)
+    _potrf_info_impl!(p.potrf, 'L', F11, info, s, sqrt(real(T)(epsc)))
     m > 0 || return nothing
     F21 = view(P, (w + 1):f, 1:w)
     _trsm_impl!(p.trsm, 'R', 'L', 'C', 'N', one(T), F11, F21)
@@ -255,23 +289,27 @@ function factorize_cholesky!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector;
     plan = N.plan
     nodes = S.schedule.group_nodes
     flag = opts.user_host_interrupt
+    # cuDSS pivot_epsilon semantics on the Cholesky path (issue #127): pivots in
+    # (0, epsc] are statically perturbed to epsc so a near-singular SPD system
+    # factors with bounded growth; nonpositive/NaN pivots still fail the front
+    epsc = resolved_pivot_epsilon(opts, real(T))
     for k in eachindex(plan.sub_first)
         _poll_interrupt(flag)
         a, b = plan.sub_first[k], plan.sub_last[k]
-        factorize_subtrees!(N, S, nzval, a, b - a + 1, plan.sub_local[k])
+        factorize_subtrees!(N, S, nzval, a, b - a + 1, plan.sub_local[k], epsc)
     end
     for k in eachindex(plan.group_first)
         _poll_interrupt(flag)
         a, b = plan.group_first[k], plan.group_last[k]
         W = plan.group_width[k]
         if W > 0
-            factorize_fronts_b!(N, S, nzval, a, b - a + 1, plan.group_maxchild[k], W)
+            factorize_fronts_b!(N, S, nzval, a, b - a + 1, plan.group_maxchild[k], W, epsc)
         else
             zero_fronts!(N, S, a, b - a + 1)
             scatter_A!(N, S, nzval, a, b - a + 1)
             extend_add!(N, S, a, b - a + 1, plan.group_maxchild[k])
             for q in a:b
-                _factor_front_c!(N, S, nodes[q], p)
+                _factor_front_c!(N, S, nodes[q], p, epsc)
             end
         end
     end
diff --git a/src/numeric/front.jl b/src/numeric/front.jl
index 672a759..bc2c91b 100644
--- a/src/numeric/front.jl
+++ b/src/numeric/front.jl
@@ -36,13 +36,17 @@ end
 # column j of the unblocked Cholesky, first half: check the pivot and apply the
 # rank-1 update of the trailing block with the scaled column (column j itself is
 # not written, so every work item reads the same unscaled values)
-@inline function _front_chol_update!(L11, st, piv, j, s, li, front_ncols, ::Val{W}, ::Val{WG}) where {W, WG}
+@inline function _front_chol_update!(L11, st, piv, j, s, li, front_ncols, epsc, ::Val{W}, ::Val{WG}) where {W, WG}
     @inbounds begin
         w = front_ncols[s]
         if j <= w && st[1] == 0
             ajj = real(L11[_packed(j, j, W)])
-            if ajj > 0
-                rd = inv(sqrt(ajj))                  # recomputed per work item: no barrier needed;
+            # pivots in [0, epsc] are statically perturbed to epsc (cuDSS pivot_epsilon
+            # semantics; epsc = 0 restores plain potrf acceptance); negative or NaN stops.
+            # Deviation from cuDSS's |pivot| rule: a NEGATIVE pivot is never perturbed --
+            # MadNLP needs the failure signal for its inertia-driven regularization
+            if ajj > epsc || (epsc > 0 && ajj >= 0)
+                rd = inv(sqrt(max(ajj, epsc)))       # recomputed per work item: no barrier needed;
                 for i in (j + li):WG:w               # strictly below the diagonal, since the diagonal
                     L11[_packed(i, j, W)] *= rd      # is replaced in the second half
                 end
@@ -56,11 +60,13 @@ end
 
 # column j, second half: rank-1 update of the trailing block with the scaled column
 # (multiply-add only; the packed-triangle enumeration replaces the guarded square)
-@inline function _front_chol_scale!(L11, st, piv, j, s, li, front_ncols, ::Val{W}, ::Val{WG}) where {W, WG}
+@inline function _front_chol_scale!(L11, st, piv, j, s, li, front_ncols, epsc, ::Val{W}, ::Val{WG}) where {W, WG}
     @inbounds begin
         w = front_ncols[s]
         if j <= w && st[1] == 0
-            li == 1 && (L11[_packed(j, j, W)] = sqrt(real(L11[_packed(j, j, W)])))
+            # max(., epsc) mirrors the update half's static perturbation (no-op when epsc = 0:
+            # the update half already stopped on any pivot this would change)
+            li == 1 && (L11[_packed(j, j, W)] = sqrt(max(real(L11[_packed(j, j, W)]), epsc)))
             r = w - j
             for q in (li - 1):WG:(r * (r + 1) ÷ 2 - 1)
                 a, b = _tri_decode(q, r)
@@ -147,7 +153,7 @@ Without `ASM` the panels and contribution blocks must already be assembled.
 """
 @kernel function front_cholesky_kernel!(factor, stack, info, nzval, amap, amap_ptr, amap_src, nodes, bm,
                                         front_ptr, front_nrows, front_ncols, cb_ptr, child_ptr, child_list,
-                                        relind_ptr, relind, maxchild, ::Val{ASM}, ::Val{W}, ::Val{NL},
+                                        relind_ptr, relind, maxchild, epsc, ::Val{ASM}, ::Val{W}, ::Val{NL},
                                         ::Val{WG}) where {ASM, W, NL, WG}
     @uniform TT = eltype(factor)
     @uniform RT = real(eltype(factor))
@@ -184,10 +190,10 @@ Without `ASM` the panels and contribution blocks must already be assembled.
     @synchronize
     for j in 1:W
         @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1]
-        _front_chol_update!(L11, st, piv, j, s, li, front_ncols, Val(W), Val(WG))
+        _front_chol_update!(L11, st, piv, j, s, li, front_ncols, epsc, Val(W), Val(WG))
         @synchronize
         @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1]
-        _front_chol_scale!(L11, st, piv, j, s, li, front_ncols, Val(W), Val(WG))
+        _front_chol_scale!(L11, st, piv, j, s, li, front_ncols, epsc, Val(W), Val(WG))
         @synchronize
     end
     @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1]
@@ -203,11 +209,11 @@ end
 
 function _launch_front_cholesky!(factor, stack, info, nzval, amap, amap_ptr, amap_src, nodes, bm::BatchMap, count,
                                  front_ptr, front_nrows, front_ncols, cb_ptr, child_ptr, child_list, relind_ptr,
-                                 relind, maxchild, asm::Val, ::Val{W}, ::Val{WG}) where {W, WG}
+                                 relind, maxchild, epsc, asm::Val, ::Val{W}, ::Val{WG}) where {W, WG}
     _ilog2(WG)
     kernel! = front_cholesky_kernel!(KernelAbstractions.get_backend(factor), WG)
     kernel!(factor, stack, info, nzval, amap, amap_ptr, amap_src, nodes, bm, front_ptr, front_nrows,
-            front_ncols, cb_ptr, child_ptr, child_list, relind_ptr, relind, Int(maxchild), asm, Val(W),
+            front_ncols, cb_ptr, child_ptr, child_list, relind_ptr, relind, Int(maxchild), epsc, asm, Val(W),
             Val(W * (W + 1) ÷ 2), Val(WG); ndrange = WG * count * bm.nact)
     return nothing
 end
@@ -230,13 +236,13 @@ Regime-B launch: assemble and factor the `count` fronts
 batch member with one [`front_cholesky_kernel!`](@ref) launch, one workgroup
 per (front, member); `info[s]` receives each front's status. Asynchronous.
 """
-function factorize_fronts_b!(N::Numeric, S::Symbolic, nzval::AbstractVector, first::Integer, count::Integer,
-                             maxchild::Integer, W::Integer)
+function factorize_fronts_b!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector, first::Integer, count::Integer,
+                             maxchild::Integer, W::Integer, epsc::Real = 0.0) where {T}
     count > 0 || return N
     _with_width_class(Int(W)) do w
         _launch_front_cholesky!(N.factor, N.stack, N.info, nzval, S.amap, S.amap_ptr, S.amap_src, S.group_nodes,
                                 batch_map(N; first), count, S.front_ptr, S.front_nrows, S.front_ncols, S.cb_ptr, S.child_ptr,
-                                S.child_list, S.relind_ptr, S.relind, maxchild, Val(true), w,
+                                S.child_list, S.relind_ptr, S.relind, maxchild, real(T)(epsc), Val(true), w,
                                 Val(FRONT_WORKGROUP))
     end
     return N
@@ -257,13 +263,14 @@ triangle, skipped when `cb_ptr[s] == 0`): `F11 = L11 L11ᴴ`, `F21 ← F21 L11
 function front_cholesky!(factor::AbstractVector, stack::AbstractVector, info::AbstractVector{Int32},
                          nodes::AbstractVector, first::Integer, count::Integer, front_ptr::AbstractVector,
                          front_nrows::AbstractVector, front_ncols::AbstractVector, cb_ptr::AbstractVector;
-                         width::Integer = REGIME_B_MAX_WIDTH, workgroup::Val = Val(FRONT_WORKGROUP))
+                         width::Integer = REGIME_B_MAX_WIDTH, workgroup::Val = Val(FRONT_WORKGROUP),
+                         epsc::Real = 0.0)
     count > 0 || return info
     _with_width_class(Int(width)) do w
         _launch_front_cholesky!(factor, stack, info, factor, nodes, nodes, nodes, nodes, single_batch(; first), count,
                                 front_ptr,
-                                front_nrows, front_ncols, cb_ptr, nodes, nodes, nodes, nodes, 0, Val(false), w,
-                                workgroup)
+                                front_nrows, front_ncols, cb_ptr, nodes, nodes, nodes, nodes, 0,
+                                real(eltype(factor))(epsc), Val(false), w, workgroup)
     end
     return info
 end
diff --git a/src/numeric/subtree.jl b/src/numeric/subtree.jl
index 022425d..ec82d33 100644
--- a/src/numeric/subtree.jl
+++ b/src/numeric/subtree.jl
@@ -142,14 +142,18 @@ end
 # (every work item recomputes `rd` from local memory: cheaper than sharing it through
 # a barrier, and it keeps the second half free of divisions; the diagonal itself is
 # written in the second half, so no work item reads a value another is overwriting)
-@inline function _subtree_chol_update!(buf, ctl, piv, j, li, ::Val{WG}) where {WG}
+@inline function _subtree_chol_update!(buf, ctl, piv, j, li, epsc, ::Val{WG}) where {WG}
     @inbounds begin
         if ctl[_ST_STATUS] == 0
             f = ctl[_ST_F]
             lf = ctl[_ST_LF] - 1
             ajj = real(buf[lf + _packed(j, j, f)])
-            if ajj > 0
-                rd = inv(sqrt(ajj))
+            # pivots in [0, epsc] are statically perturbed to epsc (cuDSS pivot_epsilon
+            # semantics; epsc = 0 restores plain potrf acceptance); negative or NaN stops.
+            # Deviation from cuDSS's |pivot| rule: a NEGATIVE pivot is never perturbed --
+            # MadNLP needs the failure signal for its inertia-driven regularization
+            if ajj > epsc || (epsc > 0 && ajj >= 0)
+                rd = inv(sqrt(max(ajj, epsc)))
                 for i in (j + li):WG:f
                     buf[lf + _packed(i, j, f)] *= rd
                 end
@@ -165,12 +169,13 @@ end
 # (multiply-add only: the former divide-in-the-update doubled the flop cost in f64
 # divisions; the packed-triangle enumeration halves the work items of the old
 # guarded full-square loop)
-@inline function _subtree_chol_scale!(buf, ctl, piv, j, li, ::Val{WG}) where {WG}
+@inline function _subtree_chol_scale!(buf, ctl, piv, j, li, epsc, ::Val{WG}) where {WG}
     @inbounds begin
         if ctl[_ST_STATUS] == 0
             f = ctl[_ST_F]
             lf = ctl[_ST_LF] - 1
-            li == 1 && (buf[lf + _packed(j, j, f)] = sqrt(real(buf[lf + _packed(j, j, f)])))
+            # max(., epsc) mirrors the update half's static perturbation (no-op when epsc = 0)
+            li == 1 && (buf[lf + _packed(j, j, f)] = sqrt(max(real(buf[lf + _packed(j, j, f)]), epsc)))
             r = f - j
             for q in (li - 1):WG:(r * (r + 1) ÷ 2 - 1)
                 a, b = _tri_decode(q, r)
@@ -247,7 +252,7 @@ write it to the update stack (packed lower triangle).
 @kernel function subtree_cholesky_kernel!(factor, stack, info, nzval, amap, amap_ptr, amap_src, trees, bm,
                                           subtree_ptr, subtree_nodes, front_ptr, front_nrows, front_ncols, cb_ptr,
                                           local_front, local_cb, child_ptr, child_list, relind_ptr, relind,
-                                          ::Val{NE}, ::Val{WG}) where {NE, WG}
+                                          epsc, ::Val{NE}, ::Val{WG}) where {NE, WG}
     @uniform TT = eltype(factor)
     @uniform RT = real(eltype(factor))
     @uniform IT = eltype(subtree_nodes)
@@ -278,9 +283,9 @@ write it to the update stack (packed lower triangle).
             @synchronize
         end
         for j in 1:ctl[_ST_W]
-            _subtree_chol_update!(buf, ctl, piv, j, li, Val(WG))
+            _subtree_chol_update!(buf, ctl, piv, j, li, epsc, Val(WG))
             @synchronize
-            _subtree_chol_scale!(buf, ctl, piv, j, li, Val(WG))
+            _subtree_chol_scale!(buf, ctl, piv, j, li, epsc, Val(WG))
             @synchronize
         end
         mb = _bm_gmember(bm, G)
@@ -294,13 +299,13 @@ write it to the update stack (packed lower triangle).
     end
 end
 
-function _launch_subtrees!(N::Numeric{T}, S::Symbolic, nzval, first, count, ::Val{LB}, ::Val{WG}) where {T, LB, WG}
+function _launch_subtrees!(N::Numeric{T}, S::Symbolic, nzval, first, count, epsc, ::Val{LB}, ::Val{WG}) where {T, LB, WG}
     _ilog2(WG)
     bm = batch_map(N; first)
     kernel! = subtree_cholesky_kernel!(KernelAbstractions.get_backend(N.factor), WG)
     kernel!(N.factor, N.stack, N.info, nzval, S.amap, S.amap_ptr, S.amap_src, S.group_nodes, bm,
             S.subtree_ptr, S.subtree_nodes, S.front_ptr, S.front_nrows, S.front_ncols, S.cb_ptr, S.local_front,
-            S.local_cb, S.child_ptr, S.child_list, S.relind_ptr, S.relind,
+            S.local_cb, S.child_ptr, S.child_list, S.relind_ptr, S.relind, real(T)(epsc),
             Val((LB - SUBTREE_LOCAL_RESERVE) ÷ sizeof(T)), Val(WG); ndrange = WG * count * bm.nact)
     return nothing
 end
@@ -325,10 +330,10 @@ are written, and each subtree root's contribution block goes to the update
 stack. Asynchronous.
 """
 function factorize_subtrees!(N::Numeric, S::Symbolic, nzval::AbstractVector, first::Integer, count::Integer,
-                             local_bytes::Integer)
+                             local_bytes::Integer, epsc::Real = 0.0)
     count > 0 || return N
     _with_local_bytes(Int(local_bytes)) do lb
-        _launch_subtrees!(N, S, nzval, first, count, lb, Val(SUBTREE_WORKGROUP))
+        _launch_subtrees!(N, S, nzval, first, count, epsc, lb, Val(SUBTREE_WORKGROUP))
     end
     return N
 end
diff --git a/src/options.jl b/src/options.jl
index f0fcded..c9b19cf 100644
--- a/src/options.jl
+++ b/src/options.jl
@@ -112,7 +112,7 @@ defaults below. Keyword arguments are applied through [`setparam!`](@ref), so
 | `ir_tol` | `0.0` | no early exit |
 | `pivot_type` | `PIVOT_AUTO` | Bunch–Kaufman for `S`/`H`, none for `SPD`/`HPD` |
 | `pivot_threshold` | `0.01` | pivot acceptance threshold; with `pivot_pairs = "default"` also the partner threshold of the 2×2 pairs chosen at analysis (a value set after `"analysis"` changes the in-front test, not the pairs already chosen) |
-| `pivot_epsilon` | `nothing` | [`default_pivot_epsilon`](@ref)`(T)` |
+| `pivot_epsilon` | `nothing` | [`default_pivot_epsilon`](@ref)`(T)`; LDLᵀ/LDLᴴ/LU static pivot perturbation, and on the Cholesky path (issue #127) the threshold: a pivot in `[0, pivot_epsilon]` is raised to it (negative or NaN still fails the front); `0` disables |
 | `pivot_epsilon_alg` | `PIVOT_EPSILON_DEFAULT` | |
 | `max_lu_nnz` | `-1` | no limit (any negative value) |
 | `hybrid_memory_mode` | `0` | device memory only |
diff --git a/src/reference/cholesky.jl b/src/reference/cholesky.jl
index 291e41e..3f944bc 100644
--- a/src/reference/cholesky.jl
+++ b/src/reference/cholesky.jl
@@ -1,3 +1,7 @@
+# NOTE (#127): the reference Cholesky is the UNPERTURBED ground truth -- it does not
+# apply pivot_epsilon (the device default statically perturbs pivots in [0, epsilon]).
+# Tests that plant near-zero pivots pin pivot_epsilon = 0 on the device side, and the
+# panel comparisons only use generators whose pivots sit far above the defaults.
 # CPU reference multifrontal Cholesky (PLAN §7, the oracle): plain Julia on
 # host arrays with BLAS/LAPACK, processing the supernodes of a `Symbolic` in
 # postorder with the same maps the device path uses (`amap_ptr`/`amap_src` for
diff --git a/test/matrices.jl b/test/matrices.jl
index 1dae513..d25eaa7 100644
--- a/test/matrices.jl
+++ b/test/matrices.jl
@@ -239,8 +239,9 @@ Symmetric (Hermitian for complex `T`) `n × n` matrix whose row and column `j`
 are zero: an SPD pentadiagonal matrix (diagonal 6, couplings of modulus ≤ 1)
 on the other `n - 1` indices, and a structurally zero pivot at `j` (stored as
 an explicit zero when `stored_zero = true`). Any elimination order meets a zero
-pivot at column `j`, so Cholesky reports `info = j` and LDLᵀ must perturb that
-pivot; `b = A x` is a consistent right-hand side. Setting `A[j, j]` to a
+pivot at column `j`, so Cholesky reports `info = j` under `pivot_epsilon = 0`
+(the default epsilon statically perturbs a zero pivot, #127) and LDLᵀ must
+perturb that pivot; `b = A x` is a consistent right-hand side. Setting `A[j, j]` to a
 positive value makes the matrix SPD, a negative value makes Cholesky fail at `j`.
 """
 function singular_block_matrix(::Type{T}, n::Integer, j::Integer; stored_zero::Bool = false) where {T}
diff --git a/test/test_api.jl b/test/test_api.jl
index 06a230c..07bcd89 100644
--- a/test/test_api.jl
+++ b/test/test_api.jl
@@ -472,8 +472,11 @@ end
     cholesky!(solver, api_matrix(backend, tril(A), INT))
     @test !solver.fresh_factorization
     @test relres(A, to_host(solver \ bd), b) <= tol(T)
-    # check = true throws on a failed factorization; the default does not (CUDSS.jl)
-    B = random_spd(T, n, 0.03) - n * I
+    # check = true throws on a failed factorization; the default does not (CUDSS.jl).
+    # -3n I makes the FIRST diagonal pivot decisively negative (~-2n): a negative pivot
+    # is never perturbed, so the failure holds under the default pivot_epsilon too (#127)
+    # -- the old -n I shift left tiny-positive pivots that ComplexF32's 1e-5 epsilon lifted
+    B = random_spd(T, n, 0.03) - 3 * n * I
     @test thrown(() -> cholesky(api_matrix(backend, tril(B), INT); view = 'L', check = true)) isa FactorizationError
     Fb = cholesky(api_matrix(backend, tril(B), INT); view = 'L')
     @test getparam(Fb, "info") > 0
diff --git a/test/test_numeric_cholesky_a.jl b/test/test_numeric_cholesky_a.jl
index a3db07b..892e1dc 100644
--- a/test/test_numeric_cholesky_a.jl
+++ b/test/test_numeric_cholesky_a.jl
@@ -201,7 +201,10 @@ end
             s = S.partition.col2sn[S.partition.iperm[j]]
             @test S.schedule.regime[s] == SDS.REGIME_A  # the failing front is in a subtree
             @test info_ref == j
-            @test SDS.factorize!(Nd, Sd, nz) == j
+            # pivot_epsilon = 0: this testset is about FAILURE REPORTING, and the default
+            # epsilon would statically perturb the d = 0 pivot (cuDSS semantics, #127 --
+            # covered by the pivot_epsilon testset below); the reference is unperturbed
+            @test SDS.factorize!(Nd, Sd, nz; opts = Options(pivot_epsilon = 0.0)) == j
             stats = to_host(Nd.stats)
             @test stats[s * SDS.FRONT_STATS_FIELDS] == Nr.stats[s * SDS.FRONT_STATS_FIELDS] ==
                   S.partition.iperm[j] - S.partition.super_ptr[s] + 1
@@ -215,3 +218,75 @@ end
     @test info_ref == 4
     @test SDS.factorize!(Nd, Sd, nz) == 4
 end
+
+@testset "pivot_epsilon static perturbation ($(backend_name(backend)), $T)" for backend in BACKENDS, T in ELTYPES
+    # Cholesky pivot 2 of [1 1; 1 1] is exactly zero: with pivot_epsilon = eps it is
+    # statically perturbed to eps and the factorization succeeds with bounded growth
+    # (issue #127, the cuDSS semantics the LDL/LU paths already follow); with
+    # pivot_epsilon = 0 the factorization stops there, and a negative pivot always stops.
+    # "algo2" + no subtree budgets forces regime C (the vendor-potrf clamp); the default
+    # options exercise the kernel paths (regime A here: tiny subtree front).
+    # elimination makes pivot 2 of [1 1; 1 1 + d] exactly d: d = 0 (zero pivot) and
+    # d = epsilon/4 (tiny positive pivot, the near-singular-KKT mechanism of #127).
+    # The KERNEL routes statically perturb pivots in (0, epsilon] to epsilon (cuDSS
+    # semantics) and perturb/accept the zero pivot; the vendor-POTRF route cannot
+    # perturb mid-elimination, so it DETECTS instead: potrf rejects the nonpositive
+    # pivot itself, and the sqrt(epsilon) factor-diagonal floor converts a silently
+    # accepted tiny pivot into a nonzero info. Either way MadNLP gets a signal.
+    Random.seed!(666)
+    epsd = SparseDirectSolver.default_pivot_epsilon(T)
+    dt = T(epsd / 4)
+    b = T[1, 0]
+    mk(M, alg_c; eps0 = false) = begin
+        solver = DirectSolver(api_matrix(backend, tril(M), Int32), spd_structure(T), 'L')
+        alg_c && setparam!(solver, "factorization_alg", "algo2")
+        alg_c && (solver.options.subtree_budgets = Int[])            # no regime A: potrf path
+        eps0 && setparam!(solver, "pivot_epsilon", 0.0)
+        solver
+    end
+    fact(M, alg_c; eps0 = false) = begin
+        solver = mk(M, alg_c; eps0)
+        x = to_device(backend, zeros(T, 2))
+        bd = to_device(backend, copy(b))
+        execute!("analysis", solver, x, bd; asynchronous = false)
+        execute!("factorization", solver, x, bd; asynchronous = false)
+        info = getparam(solver, "info")
+        info == 0 && execute!("solve", solver, x, bd; asynchronous = false)
+        (info, to_host(x))
+    end
+    Az = sparse(T[1 1; 1 1])              # pivot 2 = 0 (emerges during elimination)
+    At = sparse(T[1 1; 1 1 + dt])         # pivot 2 = dt in (0, epsilon] (during elimination)
+    As = sparse(T[1 0; 0 dt])             # tiny pivot STORED on the input diagonal
+    An = sparse(T[1 2; 2 1])              # pivot 2 < 0
+    for (rname, alg_c, perturbs) in (("kernel regimes", false, true), ("regime C / potrf", true, false))
+        @testset "$rname" begin
+            it, xt = fact(At, alg_c)
+            if perturbs
+                @test it == 0                                        # tiny pivot perturbed
+                @test all(isfinite, xt)                              # bounded growth (system is ~1/eps conditioned)
+                iz, xz = fact(Az, alg_c)
+                @test iz == 0                                        # zero pivot perturbed
+                @test all(isfinite, xz)
+            else
+                @test it != 0                                        # tiny pivot DETECTED by the epsilon floor
+                iz, _ = fact(Az, alg_c)
+                @test iz != 0                                        # zero pivot rejected by potrf itself
+            end
+            # a tiny pivot stored on the INPUT diagonal: the kernel routes perturb it, and
+            # the potrf route's pre-clamp lifts it so potrf succeeds with the factored
+            # diagonal exactly sqrt(epsilon) -- which the floor must accept (equality edge)
+            is_, xs = fact(As, alg_c)
+            @test is_ == 0
+            @test all(isfinite, xs)
+            # pivot_epsilon = 0 restores plain acceptance: the tiny positive pivot passes
+            # silently (the pre-#127 behavior), the zero pivot fails
+            i0t, _ = fact(At, alg_c; eps0 = true)
+            @test i0t == 0
+            i0z, _ = fact(Az, alg_c; eps0 = true)
+            @test i0z != 0
+            # a negative pivot is never perturbed, with or without the epsilon
+            inn, _ = fact(An, alg_c)
+            @test inn != 0
+        end
+    end
+end
diff --git a/test/test_numeric_cholesky_b.jl b/test/test_numeric_cholesky_b.jl
index 6f2b53f..3e0b942 100644
--- a/test/test_numeric_cholesky_b.jl
+++ b/test/test_numeric_cholesky_b.jl
@@ -211,7 +211,9 @@ end
             s = S.partition.col2sn[S.partition.iperm[j]]
             @test S.schedule.regime[s] == SDS.REGIME_B  # the failing front is a fused one
             @test info_ref == j
-            @test SDS.factorize!(Nd, Sd, nz) == j
+            # pivot_epsilon = 0: failure REPORTING is under test; the default epsilon
+            # statically perturbs the d = 0 pivot (#127; see cholesky_a's epsilon testset)
+            @test SDS.factorize!(Nd, Sd, nz; opts = Options(pivot_epsilon = 0.0)) == j
             stats = to_host(Nd.stats)
             @test stats[s * SDS.FRONT_STATS_FIELDS] == Nr.stats[s * SDS.FRONT_STATS_FIELDS] ==
                   S.partition.iperm[j] - S.partition.super_ptr[s] + 1
diff --git a/test/test_numeric_cholesky_c.jl b/test/test_numeric_cholesky_c.jl
index 154723c..2781aed 100644
--- a/test/test_numeric_cholesky_c.jl
+++ b/test/test_numeric_cholesky_c.jl
@@ -136,7 +136,9 @@ end
             A[j, j] = d
             S, Nr, info_ref, Sd, Nd, nz = numeric_c_setup(backend, A; opts)
             @test info_ref == j
-            @test SDS.factorize!(Nd, Sd, nz) == j
+            # pivot_epsilon = 0: failure REPORTING is under test; the default epsilon's
+            # input-diagonal clamp would lift the stored-zero pivot before potrf (#127)
+            @test SDS.factorize!(Nd, Sd, nz; opts = Options(pivot_epsilon = 0.0)) == j
             s = S.partition.col2sn[S.partition.iperm[j]]
             stats = to_host(Nd.stats)
             @test stats[s * SDS.FRONT_STATS_FIELDS] == Nr.stats[s * SDS.FRONT_STATS_FIELDS] ==