diff --git a/src/IPM/HSD/step.jl b/src/IPM/HSD/step.jl
index c39f6b77..5b48a7bf 100644
--- a/src/IPM/HSD/step.jl
+++ b/src/IPM/HSD/step.jl
@@ -86,7 +86,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
     α = max_step_length(pt, Δ)
     γ = (one(T) - α)^2 * min(one(T) - α, params.GammaMin)
     η = one(T) - γ
-    
+
     # Mehrotra corrector
     @timeit hsd.timer "Newton" solve_newton_system!(Δ, hsd, hx, hy, h0,
         # Right-hand side of Newton system
@@ -104,8 +104,8 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
         ncor += 1
 
         # Compute extra-corrector
-        αc = compute_higher_corrector!(Δc, 
-            hsd, γ, 
+        αc = compute_higher_corrector!(Δc,
+            hsd, γ,
             hx, hy, h0,
             Δ, α_, params.CentralityOutlierThreshold
         )
@@ -121,7 +121,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
             Δ.κ   = Δc.κ
             α = αc
         end
-        
+
         if αc < T(11 // 10) * α_
             break
         end
@@ -130,7 +130,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
         #     # not enough improvement, step correcting
         #     break
         # end
-        
+
     end
 
     # Update current iterate
@@ -144,7 +144,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
     pt.τ   += α  * Δ.τ
     pt.κ   += α  * Δ.κ
     update_mu!(pt)
-    
+
     return nothing
 end
 
@@ -213,11 +213,11 @@ function solve_newton_system!(Δ::Point{T, Tv},
 
     # II. Recover Δτ, Δx, Δy
     # Compute Δτ
-    @timeit hsd.timer "ξg_" ξg_ = (ξg + ξtk / pt.τ 
+    @timeit hsd.timer "ξg_" ξg_ = (ξg + ξtk / pt.τ
         - dot((ξxzl ./ pt.xl) .* dat.lflag, dat.l .* dat.lflag)            # l'(Xl)^-1 * ξxzl
-        + dot((ξxzu ./ pt.xu) .* dat.uflag, dat.u .* dat.uflag) 
-        - dot(((pt.zl ./ pt.xl) .* ξl) .* dat.lflag, dat.l .* dat.lflag) 
-        - dot(((pt.zu ./ pt.xu) .* ξu) .* dat.uflag, dat.u .* dat.uflag)  # 
+        + dot((ξxzu ./ pt.xu) .* dat.uflag, dat.u .* dat.uflag)
+        - dot(((pt.zl ./ pt.xl) .* ξl) .* dat.lflag, dat.l .* dat.lflag)
+        - dot(((pt.zu ./ pt.xu) .* ξu) .* dat.uflag, dat.u .* dat.uflag)  #
     )
 
     @timeit hsd.timer "Δτ" Δ.τ = (
@@ -276,7 +276,7 @@ function max_step_length(x::Vector{T}, dx::Vector{T}) where{T}
     n == size(dx, 1) || throw(DimensionMismatch())
     a = T(Inf)
 
-    #=@inbounds=# for i in Base.OneTo(n)
+    @inbounds for i in 1:n
         if dx[i] < zero(T)
             if (-x[i] / dx[i]) < a
                 a = (-x[i] / dx[i])
@@ -299,7 +299,7 @@ function max_step_length(pt::Point{T, Tv}, δ::Point{T, Tv}) where{T, Tv<:Abstra
 
     at = δ.τ < zero(T) ? (-pt.τ / δ.τ) : oneunit(T)
     ak = δ.κ < zero(T) ? (-pt.κ / δ.κ) : oneunit(T)
-    
+
     α = min(one(T), axl, axu, azl, azu, at, ak)
 
     return α
@@ -316,7 +316,7 @@ Requires the solution of one Newton system.
 # Arguments
 - `Δc`: Corrected search direction, modified in-place
 - `hsd`: The HSD optimizer
-- `γ`: 
+- `γ`:
 - `hx, hy, h0`: Terms obtained from the preliminary augmented system solve
 - `Δ`: Current predictor direction
 - `α`: Maximum step length in predictor direction
@@ -398,4 +398,4 @@ function compute_higher_corrector!(Δc::Point{T, Tv},
     # Compute corrected step-length
     αc = max_step_length(pt, Δc)
     return αc
-end
\ No newline at end of file
+end
diff --git a/src/IPM/MPC/MPC.jl b/src/IPM/MPC/MPC.jl
index fc56bf1e..eec253e6 100644
--- a/src/IPM/MPC/MPC.jl
+++ b/src/IPM/MPC/MPC.jl
@@ -16,7 +16,6 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
     primal_status::SolutionStatus     # Primal solution status
     dual_status::SolutionStatus       # Dual   solution status
 
-
     primal_objective::T  # Primal bound: c'x
     dual_objective::T    # Dual bound: b'y + l' zl - u'zu
 
@@ -28,7 +27,7 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
     pt::Point{T, Tv}       # Current primal-dual iterate
     res::Residuals{T, Tv}  # Residuals at current iterate
 
-    Δ::Point{T, Tv}   # Predictor direction
+    Δ::Point{T, Tv}   # Predictor
     Δc::Point{T, Tv}  # Corrector
 
     # Step sizes
@@ -36,7 +35,14 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
     αd::T
 
     # Newton system RHS
-
+    ξp::Tv
+    ξl::Tv
+    ξu::Tv
+    ξd::Tv
+    ξxzl::Tv
+    ξxzu::Tv
+
+    # KKT solver
     kkt::Tk
     regP::Tv  # Primal regularization
     regD::Tv  # Dual regularization
@@ -44,11 +50,11 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
     function MPC(
         dat::IPMData{T, Tv, Tb, Ta}, kkt_options::KKTOptions{T}
     ) where{T, Tv<:AbstractVector{T}, Tb<:AbstractVector{Bool}, Ta<:AbstractMatrix{T}}
-        
+
         m, n = dat.nrow, dat.ncol
         p = sum(dat.lflag) + sum(dat.uflag)
 
-        # Allocate some memory
+        # Working memory
         pt  = Point{T, Tv}(m, n, p, hflag=false)
         res = Residuals(
             tzeros(Tv, m), tzeros(Tv, n), tzeros(Tv, n),
@@ -58,6 +64,14 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
         Δ  = Point{T, Tv}(m, n, p, hflag=false)
         Δc = Point{T, Tv}(m, n, p, hflag=false)
 
+        # Newton RHS
+        ξp = tzeros(Tv, m)
+        ξl = tzeros(Tv, n)
+        ξu = tzeros(Tv, n)
+        ξd = tzeros(Tv, n)
+        ξxzl = tzeros(Tv, n)
+        ξxzu = tzeros(Tv, n)
+
         # Initial regularizations
         regP = tones(Tv, n)
         regD = tones(Tv, m)
@@ -70,6 +84,7 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
             T(Inf), T(-Inf),
             TimerOutput(),
             pt, res, Δ, Δc, zero(T), zero(T),
+            ξp, ξl, ξu, ξd, ξxzl, ξxzu,
             kkt, regP, regD
         )
     end
@@ -160,7 +175,7 @@ function update_solver_status!(mpc::MPC{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher
     else
         mpc.dual_status = Sln_Unknown
     end
-    
+
     # Check for optimal solution
     if ρp <= ϵp && ρd <= ϵd && ρg <= ϵg
         mpc.primal_status = Sln_Optimal
@@ -168,7 +183,7 @@ function update_solver_status!(mpc::MPC{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher
         mpc.solver_status = Trm_Optimal
         return nothing
     end
-    
+
     # TODO: Primal/Dual infeasibility detection
     # Check for infeasibility certificates
     if max(
@@ -257,12 +272,12 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}
         if params.OutputLevel > 0
             # Display log
             @printf "%4d" mpc.niter
-            
+
             # Objectives
             ϵ = dat.objsense ? one(T) : -one(T)
             @printf "  %+14.7e" ϵ * mpc.primal_objective
             @printf "  %+14.7e" ϵ * mpc.dual_objective
-            
+
             # Residuals
             @printf "  %8.2e" max(mpc.res.rp_nrm, mpc.res.rl_nrm, mpc.res.ru_nrm)
             @printf " %8.2e" mpc.res.rd_nrm
@@ -295,14 +310,14 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}
             || mpc.solver_status == Trm_DualInfeasible
         )
             break
-        elseif mpc.niter >= params.IterationsLimit 
+        elseif mpc.niter >= params.IterationsLimit
             mpc.solver_status = Trm_IterationLimit
             break
         elseif ttot >= params.TimeLimit
             mpc.solver_status = Trm_TimeLimit
             break
         end
-        
+
 
         # TODO: step
         # For now, include the factorization in the step function
@@ -314,7 +329,7 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}
             if isa(err, PosDefException) || isa(err, SingularException)
                 # Numerical trouble while computing the factorization
                 mpc.solver_status = Trm_NumericalProblem
-    
+
             elseif isa(err, OutOfMemoryError)
                 # Out of memory
                 mpc.solver_status = Trm_MemoryLimit
@@ -328,16 +343,13 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}
 
             break
         end
-
         mpc.niter += 1
-
     end
-    
+
     # TODO: print message based on termination status
     params.OutputLevel > 0 && println("Solver exited with status $((mpc.solver_status))")
 
     return nothing
-
 end
 
 function compute_starting_point(mpc::MPC{T}) where{T}
@@ -351,7 +363,7 @@ function compute_starting_point(mpc::MPC{T}) where{T}
     # Get initial iterate
     KKT.solve!(zeros(T, n), pt.y, mpc.kkt, false .* mpc.dat.b, mpc.dat.c)  # For y
     KKT.solve!(pt.x, zeros(T, m), mpc.kkt, mpc.dat.b, false .* mpc.dat.c)  # For x
-    
+
     # I. Recover positive primal-dual coordinates
     δx = one(T) + max(
         zero(T),
@@ -375,11 +387,11 @@ function compute_starting_point(mpc::MPC{T}) where{T}
     =#
     pt.zl .= ( z ./ (dat.lflag + dat.uflag)) .* dat.lflag
     pt.zu .= (-z ./ (dat.lflag + dat.uflag)) .* dat.uflag
-    
+
     δz = one(T) + max(zero(T), (-3 // 2) * minimum(pt.zl), (-3 // 2) * minimum(pt.zu))
     pt.zl[dat.lflag] .+= δz
     pt.zu[dat.uflag] .+= δz
-    
+
     mpc.pt.τ   = one(T)
     mpc.pt.κ   = zero(T)
 
@@ -397,4 +409,4 @@ function compute_starting_point(mpc::MPC{T}) where{T}
     update_mu!(mpc.pt)
 
     return nothing
-end
\ No newline at end of file
+end
diff --git a/src/IPM/MPC/step.jl b/src/IPM/MPC/step.jl
index 9fe3c720..f1440aea 100644
--- a/src/IPM/MPC/step.jl
+++ b/src/IPM/MPC/step.jl
@@ -51,44 +51,74 @@ function compute_step!(mpc::MPC{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
     nbump < 3 || throw(PosDefException(0))  # factorization could not be saved
 
     # II. Compute search direction
-    Δ = mpc.Δ
-    
-    # Affine-scaling direction (predictor)
-    @timeit mpc.timer "Newton" solve_newton_system!(Δ, mpc,
-        # Right-hand side of Newton system
-        res.rp, res.rl, res.ru, res.rd,
-        -(pt.xl .* pt.zl) .* dat.lflag,
-        -(pt.xu .* pt.zu) .* dat.uflag,
-    )
+    Δ  = mpc.Δ
+    Δc = mpc.Δc
 
-    # Step length for affine-scaling direction
-    αp, αd = max_step_length_pd(pt, Δ)
-    μₐ = (
-        dot((pt.xl + αp .* Δ.xl) .* dat.lflag, pt.zl + αd .* Δ.zl)
-        + dot((pt.xu + αp .* Δ.xu) .* dat.uflag, pt.zu + αd .* Δ.zu)
-    ) / pt.p
-    σ = clamp((μₐ / pt.μ)^3, sqrt(eps(T)), one(T) - sqrt(eps(T)))
-    
-    # Mehrotra corrector
-    @timeit mpc.timer "Newton" solve_newton_system!(Δ, mpc,
-        # Right-hand side of Newton system
-        res.rp, res.rl, res.ru, res.rd,
-        (-(pt.xl .* pt.zl) .+ σ * pt.μ .- Δ.xl .* Δ.zl) .* dat.lflag,
-        (-(pt.xu .* pt.zu) .+ σ * pt.μ .- Δ.xu .* Δ.zu) .* dat.uflag,
-    )
-    αp, αd = max_step_length_pd(pt, Δ)
+    # Affine-scaling direction and associated step size
+    @timeit mpc.timer "Predictor" compute_predictor!(mpc::MPC)
+    mpc.αp, mpc.αd = max_step_length_pd(mpc.pt, mpc.Δ)
+
+    # TODO: if step size is large enough, skip corrector
+
+    # Corrector
+    @timeit mpc.timer "Corrector" compute_corrector!(mpc::MPC)
+    mpc.αp, mpc.αd = max_step_length_pd(mpc.pt, mpc.Δc)
+    # TODO: the following is not needed if there are no additional corrections
+    copyto!(Δ.x, Δc.x)
+    copyto!(Δ.xl, Δc.xl)
+    copyto!(Δ.xu, Δc.xu)
+    copyto!(Δ.y, Δc.y)
+    copyto!(Δ.zl, Δc.zl)
+    copyto!(Δ.zu, Δc.zu)
+
+    # Extra centrality corrections
+    ncor = 0
+    ncor_max = params.CorrectionLimit
+
+    # Zero out the Newton RHS. This only needs to be done once.
+    # TODO: not needed if no additional corrections
+    rmul!(mpc.ξp, zero(T))
+    rmul!(mpc.ξl, zero(T))
+    rmul!(mpc.ξu, zero(T))
+    rmul!(mpc.ξd, zero(T))
+
+    @timeit mpc.timer "Extra corr" while ncor < ncor_max
+        compute_extra_correction!(mpc)
+
+        # TODO: function to compute step size given Δ and Δc
+        # This would avoid copying data around
+        αp_c, αd_c = max_step_length_pd(mpc.pt, mpc.Δc)
+
+        if αp_c >= 1.01 * mpc.αp && αd_c >= 1.01 * mpc.αd
+            mpc.αp = αp_c
+            mpc.αd = αd_c
+
+            # Δ ⟵ Δc
+            copyto!(Δ.x, Δc.x)
+            copyto!(Δ.xl, Δc.xl)
+            copyto!(Δ.xu, Δc.xu)
+            copyto!(Δ.y, Δc.y)
+            copyto!(Δ.zl, Δc.zl)
+            copyto!(Δ.zu, Δc.zu)
+
+            ncor += 1
+        else
+            # Not enough improvement: abort
+            break
+        end
+    end
 
     # Update current iterate
-    αp *= params.StepDampFactor
-    αd *= params.StepDampFactor
-    pt.x  .+= αp .* Δ.x
-    pt.xl .+= αp .* Δ.xl
-    pt.xu .+= αp .* Δ.xu
-    pt.y  .+= αd .* Δ.y
-    pt.zl .+= αd .* Δ.zl
-    pt.zu .+= αd .* Δ.zu
+    mpc.αp *= params.StepDampFactor
+    mpc.αd *= params.StepDampFactor
+    pt.x  .+= mpc.αp .* Δ.x
+    pt.xl .+= mpc.αp .* Δ.xl
+    pt.xu .+= mpc.αp .* Δ.xu
+    pt.y  .+= mpc.αd .* Δ.y
+    pt.zl .+= mpc.αd .* Δ.zl
+    pt.zu .+= mpc.αd .* Δ.zu
     update_mu!(pt)
-    
+
     return nothing
 end
 
@@ -167,10 +197,10 @@ function solve_newton_system!(Δ::Point{T, Tv},
 
     # Check Newton residuals
     # @printf "Newton residuals:\n"
-    # @printf "|rp|   = %16.8e\n" norm(dat.A * Δ.x + mpc.regD .* Δ.y - ξp, Inf)
+    # @printf "|rp|   = %16.8e\n" norm(dat.A * Δ.x - ξp, Inf)
     # @printf "|rl|   = %16.8e\n" norm((Δ.x - Δ.xl) .* dat.lflag - ξl, Inf)
     # @printf "|ru|   = %16.8e\n" norm((Δ.x + Δ.xu) .* dat.uflag - ξu, Inf)
-    # @printf "|rd|   = %16.8e\n" norm(-mpc.regP .* Δ.x + dat.A'Δ.y + Δ.zl - Δ.zu - ξd, Inf)
+    # @printf "|rd|   = %16.8e\n" norm(dat.A'Δ.y + Δ.zl - Δ.zu - ξd, Inf)
     # @printf "|rxzl| = %16.8e\n" norm(pt.zl .* Δ.xl + pt.xl .* Δ.zl - ξxzl, Inf)
     # @printf "|rxzu| = %16.8e\n" norm(pt.zu .* Δ.xu + pt.xu .* Δ.zu - ξxzu, Inf)
 
@@ -192,4 +222,139 @@ function max_step_length_pd(pt::Point{T, Tv}, δ::Point{T, Tv}) where{T, Tv<:Abs
     αd = min(one(T), azl, azu)
 
     return αp, αd
-end
\ No newline at end of file
+end
+
+
+"""
+    compute_predictor!(mpc::MPC) -> Nothing
+"""
+function compute_predictor!(mpc::MPC)
+
+    # Newton RHS
+    copyto!(mpc.ξp, mpc.res.rp)
+    copyto!(mpc.ξl, mpc.res.rl)
+    copyto!(mpc.ξu, mpc.res.ru)
+    copyto!(mpc.ξd, mpc.res.rd)
+    @. mpc.ξxzl = -(mpc.pt.xl .* mpc.pt.zl) .* mpc.dat.lflag
+    @. mpc.ξxzu = -(mpc.pt.xu .* mpc.pt.zu) .* mpc.dat.uflag
+
+    # Compute affine-scaling direction
+    @timeit mpc.timer "Newton" solve_newton_system!(mpc.Δ, mpc,
+        mpc.ξp, mpc.ξl, mpc.ξu, mpc.ξd, mpc.ξxzl, mpc.ξxzu
+    )
+
+    # TODO: check Newton system residuals, perform iterative refinement if needed
+    return nothing
+end
+
+"""
+    compute_corrector!(mpc::MPC) -> Nothing
+"""
+function compute_corrector!(mpc::MPC{T, Tv}) where{T, Tv<:AbstractVector{T}}
+    dat = mpc.dat
+    pt = mpc.pt
+    Δ = mpc.Δ
+    Δc = mpc.Δc
+
+    # Step length for affine-scaling direction
+    αp_aff, αd_aff = mpc.αp, mpc.αd
+    μₐ = (
+        dot((pt.xl + αp_aff .* Δ.xl) .* dat.lflag, pt.zl + αd_aff .* Δ.zl)
+        + dot((pt.xu + αp_aff .* Δ.xu) .* dat.uflag, pt.zu + αd_aff .* Δ.zu)
+    ) / pt.p
+    σ = clamp((μₐ / pt.μ)^3, sqrt(eps(T)), one(T) - sqrt(eps(T)))
+
+    # Newton RHS
+    # compute_predictor! was called ⟹ ξp, ξl, ξu, ξd are already set
+    @. mpc.ξxzl = (σ * pt.μ .- Δ.xl .* Δ.zl .- pt.xl .* pt.zl) .* dat.lflag
+    @. mpc.ξxzu = (σ * pt.μ .- Δ.xu .* Δ.zu .- pt.xu .* pt.zu) .* dat.uflag
+
+    # Compute corrector
+    @timeit mpc.timer "Newton" solve_newton_system!(mpc.Δc, mpc,
+        mpc.ξp, mpc.ξl, mpc.ξu, mpc.ξd, mpc.ξxzl, mpc.ξxzu
+    )
+
+    # TODO: check Newton system residuals, perform iterative refinement if needed
+    return nothing
+end
+
+"""
+    compute_extra_correction!(mpc) -> Nothing
+"""
+function compute_extra_correction!(mpc::MPC{T, Tv};
+    δ::T = T(3 // 10),
+    γ::T = T(1 // 10),
+) where{T, Tv<:AbstractVector{T}}
+    pt = mpc.pt
+    Δ  = mpc.Δ
+    Δc = mpc.Δc
+    dat = mpc.dat
+
+    # Tentative step sizes and centrality parameter
+    αp, αd = mpc.αp, mpc.αd
+    αp_ = min(αp + δ, one(T))
+    αd_ = min(αd + δ, one(T))
+
+    g  = dot(pt.xl, pt.zl) + dot(pt.xu, pt.zu)
+    gₐ = dot((pt.xl + mpc.αp .* Δ.xl) .* dat.lflag, pt.zl + mpc.αd .* Δ.zl) +
+        dot((pt.xu + mpc.αp .* Δ.xu) .* dat.uflag, pt.zu + mpc.αd .* Δ.zu)
+    μ = (gₐ / g) * (gₐ / g) * (gₐ / pt.p)
+
+    # Newton RHS
+    # ξp, ξl, ξu, ξd are already at zero
+    @timeit mpc.timer "target" begin
+        compute_target!(mpc.ξxzl, pt.xl, Δ.xl, pt.zl, Δ.zl, αp_, αd_, γ, μ)
+        compute_target!(mpc.ξxzu, pt.xu, Δ.xu, pt.zu, Δ.zu, αp_, αd_, γ, μ)
+    end
+
+    @timeit mpc.timer "Newton" solve_newton_system!(Δc, mpc,
+        mpc.ξp, mpc.ξl, mpc.ξu, mpc.ξd, mpc.ξxzl, mpc.ξxzu
+    )
+
+    # Δc ⟵ Δp + Δc
+    axpy!(one(T), Δ.x, Δc.x)
+    axpy!(one(T), Δ.xl, Δc.xl)
+    axpy!(one(T), Δ.xu, Δc.xu)
+    axpy!(one(T), Δ.y, Δc.y)
+    axpy!(one(T), Δ.zl, Δc.zl)
+    axpy!(one(T), Δ.zu, Δc.zu)
+
+    # TODO: check Newton residuals
+    return nothing
+end
+
+"""
+    compute_target!(t, x, z, γ, μ)
+
+Compute centrality target.
+"""
+function compute_target!(
+    t::Vector{T},
+    x::Vector{T},
+    δx::Vector{T},
+    z::Vector{T},
+    δz::Vector{T},
+    αp::T,
+    αd::T,
+    γ::T,
+    μ::T
+) where{T}
+
+    n = length(t)
+
+    tmin = μ * γ
+    tmax = μ / γ
+
+    @inbounds for j in 1:n
+        v = (x[j] + αp * δx[j]) * (z[j] + αd * δz[j])
+        if v < tmin
+            t[j] = tmin - v
+        elseif v > tmax
+            t[j] = tmax - v
+        else
+            t[j] = zero(T)
+        end
+    end
+
+    return nothing
+end
diff --git a/src/IPM/options.jl b/src/IPM/options.jl
index 4e5636ae..8792718c 100644
--- a/src/IPM/options.jl
+++ b/src/IPM/options.jl
@@ -13,7 +13,7 @@ Base.@kwdef mutable struct IPMOptions{T}
     ToleranceIFeas::T = sqrt(eps(T))  # infeasibility
 
     # Algorithmic parameters
-    CorrectionLimit::Int = 5  # Maximum number of centrality corrections
+    CorrectionLimit::Int = 3  # Maximum number of centrality corrections
     StepDampFactor::T = T(9_995 // 10_000)  # Damp step size by this much
     GammaMin::T = T(1 // 10)
     CentralityOutlierThreshold::T = T(1 // 10)  # Relative threshold for centrality outliers
@@ -22,4 +22,4 @@ Base.@kwdef mutable struct IPMOptions{T}
     DRegMin::T = sqrt(eps(T))  # dual
 
     Factory::Factory{<:AbstractIPMOptimizer} = Factory(HSD)
-end
\ No newline at end of file
+end