diff --git a/src/qdldl/qdldl.rs b/src/qdldl/qdldl.rs
index b36b82e6..476f2500 100644
--- a/src/qdldl/qdldl.rs
+++ b/src/qdldl/qdldl.rs
@@ -85,6 +85,168 @@ pub struct QDLDLFactorisation<T = f64> {
     workspace: QDLDLWorkspace<T>,
     /// true if factorisation is symbolic only
     is_symbolic: bool,
+    /// workspace for `solve_refined`, allocated on first use
+    ir_work: Option<RefinementWorkspace<T>>,
+    /// both-triangles copy of the internal matrix, for refinement
+    /// residuals; built on first use and refreshed when values change
+    sym: Option<SymmetricCopy<T>>,
+    /// true when `sym`'s values are stale w.r.t. the internal matrix
+    sym_stale: bool,
+}
+
+// working vectors for solve_refined, all in permuted coordinates
+#[derive(Debug)]
+struct RefinementWorkspace<T> {
+    x: Vec<T>,  // current solution
+    dx: Vec<T>, // refinement step / trial solution
+    e: Vec<T>,  // residual
+}
+
+// A full (both triangles) copy of the internally permuted matrix.
+//
+// Refinement residuals `e = b - Ax` were computed from the upper triangle
+// alone, which requires exploiting symmetry with two scattered
+// read-modify-writes per off-diagonal nonzero (into `y[row]` and
+// `y[col]`).   Holding both triangles instead costs about twice the memory
+// but removes the scatter entirely: because A is symmetric, its CSC arrays
+// read as CSR describe the same matrix, so row i of A is exactly what is
+// stored as "column i".   The residual is then a sequence of independent
+// sparse dot products -- one contiguous pass over the values, gathered
+// reads of x, and an accumulator in a register.
+//
+// `sym_to_triu` gives, for each nonzero here, the index of the
+// upper-triangular source entry it takes its value from: an off-diagonal
+// source is referenced twice (once per triangle), a diagonal source once.
+// Mapping this direction rather than source-to-destination makes the value
+// refresh a gather with sequential stores instead of a scatter.
+#[derive(Debug)]
+struct SymmetricCopy<T> {
+    A: CscMatrix<T>,
+    sym_to_triu: Vec<u32>,
+}
+
+// Builds the both-triangles copy and the value map from an upper
+// triangular source.   Returns None if the pattern is too large for the
+// u32 index map, in which case refinement falls back to the triangular
+// symmetric matrix-vector product.
+fn _build_symmetric_copy<T: FloatT>(triu: &CscMatrix<T>) -> Option<SymmetricCopy<T>> {
+    let n = triu.ncols();
+    let nnz_triu = triu.nzval.len();
+    let mut ndiag = 0;
+    for j in 0..n {
+        for &i in &triu.rowval[triu.colptr[j]..triu.colptr[j + 1]] {
+            if i == j {
+                ndiag += 1;
+            }
+        }
+    }
+    let nnz_sym = 2 * nnz_triu - ndiag;
+    if nnz_sym >= u32::MAX as usize {
+        return None;
+    }
+
+    // column counts: entry (i,j), i <= j, lands in column j and, when
+    // off-diagonal, also in column i
+    let mut colptr = vec![0usize; n + 1];
+    for j in 0..n {
+        for &i in &triu.rowval[triu.colptr[j]..triu.colptr[j + 1]] {
+            colptr[j + 1] += 1;
+            if i != j {
+                colptr[i + 1] += 1;
+            }
+        }
+    }
+    for j in 0..n {
+        colptr[j + 1] += colptr[j];
+    }
+
+    let mut rowval = vec![0usize; nnz_sym];
+    let mut next = colptr[0..n].to_vec();
+    let mut sym_to_triu = vec![0u32; nnz_sym];
+    for j in 0..n {
+        let base = triu.colptr[j];
+        for (t, &i) in triu.rowval[base..triu.colptr[j + 1]].iter().enumerate() {
+            let k = base + t;
+            let pj = next[j];
+            rowval[pj] = i;
+            sym_to_triu[pj] = k as u32;
+            next[j] += 1;
+            if i != j {
+                let pi = next[i];
+                rowval[pi] = j;
+                sym_to_triu[pi] = k as u32;
+                next[i] += 1;
+            }
+        }
+    }
+
+    let A = CscMatrix {
+        m: n,
+        n,
+        colptr,
+        rowval,
+        nzval: vec![T::zero(); nnz_sym],
+    };
+    Some(SymmetricCopy { A, sym_to_triu })
+}
+
+// Computes e = b - Ax for symmetric A held in both-triangles form,
+// reading its CSC arrays as CSR (valid because A = Aᵀ), and returns
+// ||e||_∞.   Each row is an independent sparse dot product.
+//
+// The loop is memory bound -- one streamed value and one gathered x per
+// nonzero -- so the accumulators stay in registers.  Four are used and
+// reduced pairwise.   A single accumulator serializes each row on the
+// latency of one dependent add, and independent partial sums also carry a
+// tighter error bound than sequential summation, growing like n/k + k for
+// k accumulators rather than n (Higham, "Accuracy and Stability of
+// Numerical Algorithms", 2nd ed., 2002, §4.2, on blocked and pairwise
+// summation).
+//
+// Four rather than two or eight, measured on the portfolio problems:
+// one accumulator is ~8% slower than four; two is ~1% *faster* than four
+// but degrades one problem from AlmostSolved to InsufficientProgress,
+// consistent with its looser summation error; eight is indistinguishable
+// from four in time and needs more code.   Four is therefore the smallest
+// count that captures both the pipelining and the accuracy.
+fn _sym_residual<T: FloatT>(sym: &SymmetricCopy<T>, e: &mut [T], b: &[T], x: &[T]) -> T {
+    let (colptr, rowval, nzval) = (&sym.A.colptr, &sym.A.rowval, &sym.A.nzval);
+    let mut norme = T::zero();
+    // A NaN residual must be reported, and cannot be detected by the
+    // running maximum alone: IEEE maxNum returns the non-NaN operand, so a
+    // NaN would leave `norme` finite and let a non-finite solution be
+    // accepted as converged.   Tracked separately and folded in at the end.
+    let mut any_nan = false;
+    for (i, ei) in e.iter_mut().enumerate() {
+        let (f, l) = (colptr[i], colptr[i + 1]);
+        let (vals, cols) = (&nzval[f..l], &rowval[f..l]);
+        let nchunk = vals.len() / 4;
+
+        let mut s = [T::zero(); 4];
+        unsafe {
+            // Safety: rowval entries are column indices of a matrix with
+            // the same dimension as x, as built by _build_symmetric_copy.
+            for c in 0..nchunk {
+                let (v, k) = (&vals[4 * c..4 * c + 4], &cols[4 * c..4 * c + 4]);
+                s[0] += v[0] * *x.get_unchecked(k[0]);
+                s[1] += v[1] * *x.get_unchecked(k[1]);
+                s[2] += v[2] * *x.get_unchecked(k[2]);
+                s[3] += v[3] * *x.get_unchecked(k[3]);
+            }
+            for t in 4 * nchunk..vals.len() {
+                s[t & 3] += vals[t] * *x.get_unchecked(cols[t]);
+            }
+        }
+
+        let ri = b[i] - ((s[0] + s[1]) + (s[2] + s[3]));
+        *ei = ri;
+        any_nan |= ri.is_nan();
+        norme = T::max(norme, T::abs(ri));
+    }
+    if any_nan {
+        return T::nan();
+    }
+    norme
 }
 
 impl<T> QDLDLFactorisation<T>
@@ -137,9 +299,131 @@ where
         ipermute(b, tmp, &self.perm);
     }
 
+    /// Solves Ax = b like [`solve`](crate::qdldl::QDLDLFactorisation::solve),
+    /// then iteratively refines x against the matrix currently held in the
+    /// internal workspace (i.e. the values set through
+    /// [`update_values`](crate::qdldl::QDLDLFactorisation::update_values) and
+    /// friends, which may deliberately differ from the values that were
+    /// factored, e.g. by a static regularization shift).
+    ///
+    /// The refinement loop runs entirely in the internally permuted
+    /// coordinates: the permutation is applied once to `b` and once to the
+    /// returned `x`, rather than once per backsolve as with repeated calls
+    /// to [`solve`](crate::qdldl::QDLDLFactorisation::solve).  Stopping
+    /// rules: refinement ends when the residual satisfies
+    /// `norm(b - Ax) <= abstol + reltol * norm(b)` (∞-norms), when
+    /// `max_iter` passes have been made, or when a pass fails to improve
+    /// the residual norm by at least `stop_ratio` (an improving final pass
+    /// is still accepted).  Returns false if the solution or residual
+    /// became non-finite.
+    pub fn solve_refined(
+        &mut self,
+        x: &mut [T],
+        b: &[T],
+        reltol: T,
+        abstol: T,
+        max_iter: u32,
+        stop_ratio: T,
+    ) -> bool {
+        assert!(!self.is_symbolic);
+        assert_eq!(b.len(), self.D.len());
+        assert_eq!(x.len(), self.D.len());
+
+        let n = self.D.len();
+
+        // build the both-triangles copy on first use, and refresh its
+        // values whenever the internal matrix has been modified since
+        if self.sym.is_none() {
+            self.sym = _build_symmetric_copy(&self.workspace.triuA);
+            self.sym_stale = true;
+        }
+        if self.sym_stale {
+            if let Some(sym) = &mut self.sym {
+                let src = &self.workspace.triuA.nzval;
+                for (dst, &k) in zip(&mut sym.A.nzval, &sym.sym_to_triu) {
+                    *dst = src[k as usize];
+                }
+            }
+            self.sym_stale = false;
+        }
+
+        let work = self.ir_work.get_or_insert_with(|| RefinementWorkspace {
+            x: vec![T::zero(); n],
+            dx: vec![T::zero(); n],
+            e: vec![T::zero(); n],
+        });
+
+        // permute b once; all work below is in permuted coordinates
+        let bp = &mut self.workspace.fwork;
+        permute(bp, b, &self.perm);
+
+        let (Lp, Li, Lx) = (&self.L.colptr, &self.L.rowval, &self.L.nzval);
+        let Asym = self.workspace.triuA.sym_up();
+        let normb = bp.norm_inf();
+
+        // initial solve
+        let xp = &mut work.x;
+        xp.copy_from(bp);
+        _solve(Lp, Li, Lx, &self.Dinv, xp);
+
+        // computes e = bp - A*ξ and returns its norm.  Prefers the
+        // both-triangles copy (a pure gather; see `SymmetricCopy`), and
+        // falls back to the triangular symv if that copy is unavailable.
+        let sym = self.sym.as_ref();
+        let refine_error = |e: &mut [T], ξ: &[T]| -> T {
+            match sym {
+                Some(sym) => _sym_residual(sym, e, bp, ξ),
+                None => {
+                    e.copy_from(bp);
+                    Asym.symv(e, ξ, -T::one(), T::one());
+                    e.norm_inf()
+                }
+            }
+        };
+
+        let (e, dx) = (&mut work.e, &mut work.dx);
+        let mut norme = refine_error(e, xp);
+        if !norme.is_finite() {
+            return false;
+        }
+
+        for _ in 0..max_iter {
+            if norme <= (abstol + reltol * normb) {
+                // within tolerance.  Exit
+                break;
+            }
+            let lastnorme = norme;
+
+            // make a refinement: dx = A⁻¹e, prospective solution xp + dx
+            dx.copy_from(e);
+            _solve(Lp, Li, Lx, &self.Dinv, dx);
+            dx.axpby(T::one(), xp, T::one());
+
+            norme = refine_error(e, dx);
+            if !norme.is_finite() {
+                return false;
+            }
+
+            let improved_ratio = lastnorme / norme;
+            if improved_ratio < stop_ratio {
+                // insufficient improvement.  Exit
+                if improved_ratio > T::one() {
+                    std::mem::swap(xp, dx);
+                }
+                break;
+            }
+            std::mem::swap(xp, dx);
+        }
+
+        // undo the permutation into the output
+        ipermute(x, xp, &self.perm);
+        true
+    }
+
     /// Update a subset of the values of the matrix to be (re)factored.  See [`refactor`](crate::qdldl::QDLDLFactorisation::refactor)
     ///
     pub fn update_values(&mut self, indices: &[usize], values: &[T]) {
+        self.sym_stale = true;
         let nzval = &mut self.workspace.triuA.nzval; // post perm internal data
         let AtoPAPt = &self.workspace.AtoPAPt; //mapping from input matrix entries to triuA
 
@@ -151,6 +435,7 @@ where
     /// Update a subset of the values of the matrix to be (re)factored.  See [`refactor`](crate::qdldl::QDLDLFactorisation::refactor)
     ///
     pub fn scale_values(&mut self, indices: &[usize], scale: T) {
+        self.sym_stale = true;
         let nzval = &mut self.workspace.triuA.nzval; // post perm internal data
         let AtoPAPt = &self.workspace.AtoPAPt; //mapping from input matrix entries to triuA
 
@@ -164,6 +449,7 @@ where
     /// indicating the direction of shifts. See [`refactor`](crate::qdldl::QDLDLFactorisation::refactor)
     ///
     pub fn offset_values(&mut self, indices: &[usize], offset: T, signs: &[i8]) {
+        self.sym_stale = true;
         assert_eq!(indices.len(), signs.len());
 
         let nzval = &mut self.workspace.triuA.nzval; // post perm internal data
@@ -291,6 +577,9 @@ fn _qdldl_new<T: FloatT>(
         Dinv,
         workspace,
         is_symbolic: opts.logical,
+        ir_work: None,
+        sym: None,
+        sym_stale: true,
     })
 }
 
diff --git a/src/qdldl/test.rs b/src/qdldl/test.rs
index 020424f1..b8877d15 100644
--- a/src/qdldl/test.rs
+++ b/src/qdldl/test.rs
@@ -20,6 +20,10 @@ fn test_matrix_4x4() -> CscMatrix<f64> {
     }
 }
 
+fn inf_norm<T: FloatT>(a: &[T]) -> T {
+    a.iter().fold(T::zero(), |acc, x| T::max(acc, T::abs(*x)))
+}
+
 fn inf_norm_diff<T: FloatT>(a: &[T], b: &[T]) -> T {
     zip(a, b).fold(T::zero(), |acc, (x, y)| T::max(acc, T::abs(*x - *y)))
 }
@@ -229,6 +233,188 @@ fn test_solve_basic() {
     assert!(inf_norm_diff(&x, &b) <= 1e-8);
 }
 
+#[test]
+fn test_solve_refined() {
+    let A = test_matrix_4x4();
+    let mut factors = QDLDLFactorisation::new(&A, None).unwrap();
+    let x_true = [1., -2., 3., -4.];
+    let b = [20.0, -22.0, 32.0, -7.0];
+
+    // agrees with the plain solve
+    let mut x = [0.0; 4];
+    assert!(factors.solve_refined(&mut x, &b, 1e-13, 1e-12, 10, 5.0));
+    assert!(inf_norm_diff(&x_true, &x) <= 1e-10);
+
+    // refines against the values currently in the internal workspace,
+    // even where they differ from the factored ones: scale the matrix
+    // by 1.1 without refactoring, and refinement (contraction rate
+    // ||I - A⁻¹(1.1A)|| = 0.1 per pass) must converge to the solution
+    // of the *scaled* system using the stale factors
+    let indices: Vec<usize> = (0..A.nzval.len()).collect();
+    factors.scale_values(&indices, 1.1);
+    let mut x = [0.0; 4];
+    assert!(factors.solve_refined(&mut x, &b, 1e-13, 1e-12, 20, 5.0));
+    let x_scaled: Vec<f64> = x_true.iter().map(|v| v / 1.1).collect();
+    assert!(inf_norm_diff(&x_scaled, &x) <= 1e-10);
+}
+
+// A small quasidefinite KKT-like matrix [[diag(p) B'; B -diag(r)]] in
+// upper triangular CSC form, with the D signs it should factor with.
+#[cfg(test)]
+fn test_matrix_kkt_like(nx: usize, nz: usize, seed: u64) -> (CscMatrix<f64>, Vec<i8>) {
+    let mut state = seed;
+    let mut next = move || {
+        state = state
+            .wrapping_mul(6364136223846793005)
+            .wrapping_add(1442695040888963407);
+        ((state >> 33) as f64) / ((1u64 << 31) as f64) - 1.0
+    };
+
+    let n = nx + nz;
+    let mut cols: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
+    for j in 0..nx {
+        cols[j].push((j, 1.0 + next().abs()));
+    }
+    for j in 0..nz {
+        let col = nx + j;
+        for t in 0..3 {
+            let i = ((next().abs() * nx as f64) as usize + t * 7) % nx;
+            cols[col].push((i, next()));
+        }
+        cols[col].sort_by_key(|e| e.0);
+        cols[col].dedup_by_key(|e| e.0);
+        cols[col].push((col, -(1.0 + next().abs())));
+    }
+
+    let mut colptr = vec![0usize];
+    let (mut rowval, mut nzval) = (Vec::new(), Vec::new());
+    for c in &cols {
+        for &(i, v) in c {
+            rowval.push(i);
+            nzval.push(v);
+        }
+        colptr.push(rowval.len());
+    }
+    let mut signs = vec![1i8; n];
+    signs[nx..].fill(-1);
+    (
+        CscMatrix {
+            m: n,
+            n,
+            colptr,
+            rowval,
+            nzval,
+        },
+        signs,
+    )
+}
+
+// With exactly-representable values every product and sum below is
+// exact, so the both-triangles residual must agree with the triangular
+// symv *bitwise*.  Any disagreement would be a structural error in the
+// copy or its value map rather than a rounding difference.
+#[test]
+fn test_symmetric_copy_exact_on_integer_data() {
+    // upper triangular, small integers, some zero-free columns
+    //   [ 2  -1   0   3 ]
+    //   [    -4   5   0 ]
+    //   [         6  -2 ]
+    //   [             8 ]
+    let triu = CscMatrix {
+        m: 4,
+        n: 4,
+        colptr: vec![0, 1, 3, 5, 8],
+        rowval: vec![0, 0, 1, 1, 2, 0, 2, 3],
+        nzval: vec![2., -1., -4., 5., 6., 3., -2., 8.],
+    };
+
+    let mut sym = _build_symmetric_copy(&triu).unwrap();
+    for (dst, &k) in zip(&mut sym.A.nzval, &sym.sym_to_triu) {
+        *dst = triu.nzval[k as usize];
+    }
+    // 4 diagonal entries, 4 off-diagonal mirrored
+    assert_eq!(sym.A.nzval.len(), 2 * 8 - 4);
+
+    let x = [1., -2., 4., 8.];
+    let b = [100., -100., 50., 25.];
+
+    let mut e_ref = b;
+    triu.sym_up().symv(&mut e_ref, &x, -1.0, 1.0);
+
+    let mut e_new = [0.0; 4];
+    let norme = _sym_residual(&sym, &mut e_new, &b, &x);
+
+    assert_eq!(e_ref, e_new); // bitwise on exact data
+    assert_eq!(norme, inf_norm(&e_new));
+}
+
+// A non-finite solution must be reported rather than masked.  The running
+// maximum cannot see a NaN on its own, because IEEE maxNum returns the
+// non-NaN operand, so without explicit tracking a NaN residual would leave
+// the norm finite and the caller would accept a NaN solution as converged.
+#[test]
+fn test_sym_residual_reports_nan() {
+    let triu = CscMatrix {
+        m: 3,
+        n: 3,
+        colptr: vec![0, 1, 3, 6],
+        rowval: vec![0, 0, 1, 0, 1, 2],
+        nzval: vec![2.0, 1.0, 3.0, -1.0, 0.5, 4.0],
+    };
+    let mut sym = _build_symmetric_copy(&triu).unwrap();
+    for (dst, &k) in zip(&mut sym.A.nzval, &sym.sym_to_triu) {
+        *dst = triu.nzval[k as usize];
+    }
+    let b = [1.0, 2.0, 3.0];
+    let mut e = [0.0; 3];
+
+    // finite input: ordinary infinity norm
+    let n_ok: f64 = _sym_residual(&sym, &mut e, &b, &[1.0, 1.0, 1.0]);
+    assert!(n_ok.is_finite());
+
+    // one NaN in x taints its rows, and must be reported
+    let norme: f64 = _sym_residual(&sym, &mut e, &b, &[1.0, f64::NAN, 1.0]);
+    assert!(norme.is_nan(), "NaN residual must not be masked");
+}
+
+// The both-triangles copy must reproduce the matrix the triangular
+// symv represents, so that refinement residuals are unchanged in value
+// (they differ only in summation order).
+#[test]
+fn test_symmetric_copy_matches_triangular_symv() {
+    let (A, signs) = test_matrix_kkt_like(50, 35, 777);
+    let opts = QDLDLSettingsBuilder::<f64>::default()
+        .Dsigns(signs)
+        .build()
+        .unwrap();
+    let factors = QDLDLFactorisation::new(&A, Some(opts)).unwrap();
+
+    // the permuted internal matrix, and its both-triangles copy
+    let triu = &factors.workspace.triuA;
+    let mut sym = _build_symmetric_copy(triu).unwrap();
+    for (dst, &k) in zip(&mut sym.A.nzval, &sym.sym_to_triu) {
+        *dst = triu.nzval[k as usize];
+    }
+
+    let n = triu.ncols();
+    // every off-diagonal entry must appear on both sides
+    assert_eq!(sym.A.nzval.len(), 2 * triu.nzval.len() - n);
+
+    let x: Vec<f64> = (0..n).map(|i| ((i * 13) % 7) as f64 - 3.0).collect();
+    let b: Vec<f64> = (0..n).map(|i| ((i * 5) % 11) as f64 - 5.0).collect();
+
+    // reference: e = b - A*x via the triangular symv
+    let mut e_ref = b.clone();
+    triu.sym_up().symv(&mut e_ref, &x, -1.0, 1.0);
+
+    let mut e_new = vec![0.0; n];
+    let norme = _sym_residual(&sym, &mut e_new, &b, &x);
+
+    // same values up to summation order, and the returned norm agrees
+    assert!(inf_norm_diff(&e_ref, &e_new) <= 1e-10 * inf_norm(&e_ref).max(1.0));
+    assert!((norme - inf_norm(&e_new)).abs() <= 1e-12 * norme.max(1.0));
+}
+
 #[test]
 #[should_panic]
 fn test_solve_logical() {
diff --git a/src/solver/core/kktsolvers/direct/quasidef/directldlkktsolver.rs b/src/solver/core/kktsolvers/direct/quasidef/directldlkktsolver.rs
index bf44bc1c..1314f945 100644
--- a/src/solver/core/kktsolvers/direct/quasidef/directldlkktsolver.rs
+++ b/src/solver/core/kktsolvers/direct/quasidef/directldlkktsolver.rs
@@ -171,12 +171,28 @@ where
         lhsz: Option<&mut [T]>,
         settings: &CoreSettings<T>,
     ) -> bool {
-        self.ldlsolver.solve(&self.KKT, &mut self.x, &mut self.b);
-
         let is_success = {
             if settings.iterative_refinement_enable {
-                self.iterative_refinement(settings)
+                // backends may implement refinement internally (in their
+                // own permuted coordinates); otherwise refine here via
+                // repeated solves against our unpermuted KKT copy
+                let refined = self.ldlsolver.solve_refined(
+                    &mut self.x,
+                    &self.b,
+                    settings.iterative_refinement_reltol,
+                    settings.iterative_refinement_abstol,
+                    settings.iterative_refinement_max_iter,
+                    settings.iterative_refinement_stop_ratio,
+                );
+                match refined {
+                    Some(is_success) => is_success,
+                    None => {
+                        self.ldlsolver.solve(&self.KKT, &mut self.x, &mut self.b);
+                        self.iterative_refinement(settings)
+                    }
+                }
             } else {
+                self.ldlsolver.solve(&self.KKT, &mut self.x, &mut self.b);
                 self.x.is_finite()
             }
         };
@@ -253,11 +269,14 @@ where
         let is_success = self.ldlsolver.refactor(KKT);
 
         if settings.static_regularization_enable {
-            // put our internal copy of the KKT matrix back the way
-            // it was. Not necessary to fix the ldlsolver copy because
-            // this is only needed for our post-factorization IR scheme
-
-            _update_values_KKT(KKT, &map.diag_full, diag_kkt);
+            // put our internal copy of the KKT matrix back the way it
+            // was, and likewise any internal copy held by the ldlsolver.
+            // Both are only used post-factorization, for computing
+            // iterative refinement residuals against the unregularized
+            // matrix; the shift is recomputed and reapplied from the
+            // restored diagonal before the next refactorization.
+
+            _update_values(&mut self.ldlsolver, KKT, &map.diag_full, diag_kkt);
         }
 
         is_success
diff --git a/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/qdldl.rs b/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/qdldl.rs
index ee91a171..6b45e612 100644
--- a/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/qdldl.rs
+++ b/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/qdldl.rs
@@ -104,4 +104,23 @@ where
         self.factors.refactor().unwrap();
         self.factors.Dinv.is_finite()
     }
+
+    fn solve_refined(
+        &mut self,
+        x: &mut [T],
+        b: &[T],
+        reltol: T,
+        abstol: T,
+        max_iter: u32,
+        stop_ratio: T,
+    ) -> Option<bool> {
+        // QDLDL's internal matrix copy holds the values most recently
+        // written through update/scale/offset -- i.e. the unregularized
+        // diagonal, restored there after each refactorization -- so it can
+        // refine against the true KKT matrix in permuted coordinates.
+        Some(
+            self.factors
+                .solve_refined(x, b, reltol, abstol, max_iter, stop_ratio),
+        )
+    }
 }
diff --git a/src/solver/core/kktsolvers/direct/quasidef/mod.rs b/src/solver/core/kktsolvers/direct/quasidef/mod.rs
index 891c2418..a588aa47 100644
--- a/src/solver/core/kktsolvers/direct/quasidef/mod.rs
+++ b/src/solver/core/kktsolvers/direct/quasidef/mod.rs
@@ -23,4 +23,24 @@ pub trait DirectLDLSolver<T: FloatT>: DirectLDLSolverReqs + HasLinearSolverInfo
     fn offset_values(&mut self, index: &[usize], offset: T, signs: &[i8]);
     fn solve(&mut self, kkt: &CscMatrix<T>, x: &mut [T], b: &mut [T]);
     fn refactor(&mut self, kkt: &CscMatrix<T>) -> bool;
+
+    /// Solve with iterative refinement performed inside the backend,
+    /// against the backend's internal (unregularized) matrix copy.
+    /// Backends that maintain an internally permuted copy can run the
+    /// whole refinement loop in permuted coordinates, avoiding the
+    /// per-backsolve permutations of repeated `solve` calls.   Returns
+    /// None if the backend does not support this, in which case the
+    /// caller performs its own refinement via repeated `solve` calls.
+    #[allow(clippy::too_many_arguments)]
+    fn solve_refined(
+        &mut self,
+        _x: &mut [T],
+        _b: &[T],
+        _reltol: T,
+        _abstol: T,
+        _max_iter: u32,
+        _stop_ratio: T,
+    ) -> Option<bool> {
+        None
+    }
 }