diff --git a/ext/HSL/ma57_workspace.jl b/ext/HSL/ma57_workspace.jl
index b4421df9..e3bef54d 100644
--- a/ext/HSL/ma57_workspace.jl
+++ b/ext/HSL/ma57_workspace.jl
@@ -2,6 +2,7 @@ mutable struct PenaltyMA57Workspace{
   WP<:Ma57,
   K2<:AbstractMatrix,
   V<:AbstractVector,
+  VI<:Union{Nothing,AbstractVector},
   T<:Real,
 } <: AbstractHSLWorkspace
   M::WP
@@ -9,7 +10,9 @@ mutable struct PenaltyMA57Workspace{
   x::V
   work::V
   _qn_work::V
+  _info::Base.RefValue{BlasInt}  # For CompactBFGS LU factorization
   dx::V
+  _ipiv::VI # For CompactBFGS LU factorization
   σ::T
   n::Int
   m::Int
@@ -42,6 +45,8 @@ function construct_ma57_workspace(
     similar(u1, 4*length(u1)),
     similar(u1, 0),
     similar(u1),
+    nothing,
+    Ref{BlasInt}(0),
     zero(T),
     n,
     m,
@@ -72,6 +77,8 @@ function construct_ma57_workspace(
     similar(u1, 4*length(u1)),
     similar(u1, 2*length(u1)*H.B._mem),
     similar(u1),
+    Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem),
+    Ref{BlasInt}(0),
     zero(T),
     n,
     m,
@@ -214,12 +221,12 @@ function solve_system!(
   B = workspace.H.B
   n, m = workspace.n, workspace.m
   p = min(B._insert - 1, B._mem)
-  x1, x2, x3, y1, y2 = H.x1, H.x2, H.x3, H.y1, H.y2
+  x1, x2, x3, y1 = H.x1, H.x2, H.x3, H.y1
   Z1, Z2 = H.Z1, H.Z2
   Uk = @view B.Uk[:, 1:p]
   Vk = @view B.Vk[:, 1:p]
 
-  # Step 0: Write (#TODO: we can use easily use QRMumps instead of LDLFactorization here...)
+  # Step 0: Write
   # [B  Aᵀ] = [σI+ξI  Aᵀ] + [-U V]([U V])ᵀ
   # [A -αI] = [A     -αI] + [ 0 0]([0 0])
   # Hence,
@@ -273,79 +280,104 @@ function solve_system!(
     return
   end
 
-  # Step 3: Compute
-  # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)]
-  # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)]
-  @views mul!(y1[1:p], Uk', x1[1:n])
-  @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n])
+  if p > 0
 
+    # Step 3: Compute
+    # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)]
+    # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)]
+    @views mul!(y1[1:p], Uk', x1[1:n])
+    @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n])
 
-  # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )
-  #                                   (       [A     -αI]     )
-  # Step 4.1: Compute 
-  # Z₁ = [σI+ξI  Aᵀ]⁻¹ E = [σI+ξI  Aᵀ]⁻¹[-U V]
-  # Z₁ = [A     -αI]   E = [A     -αI]  [ 0 0]
-  Z1 .= 0
 
-  @views Z1[1:n, 1:p] .= Uk .* (-1)
-  @views Z1[1:n, (p+1):(2*p)] .= Vk
-  try
-    ma57_solve!(workspace.M, Z1, workspace._qn_work)
-  catch e
-    !(e isa HSL.Ma57Exception) && rethrow(e)
-  end
-  if any(isnan, Z1) || workspace.M.info.info[1] < 0
-    workspace.status = :failed
-    return
-  end
+    # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )
+    #                                   (       [A     -αI]     )
+    # Step 4.1: Compute 
+    # Z₁ = [σI+ξI  Aᵀ]⁻¹ E = [σI+ξI  Aᵀ]⁻¹[-U V]
+    # Z₁ = [A     -αI]   E = [A     -αI]  [ 0 0]
+    Z1 .= 0
 
-  # Step 4.2: Compute 
-  # Z₂ = FᵀZ₁ = UᵀZ₁[1:n]
-  # Z₂ = FᵀZ₁ = VᵀZ₁[1:n]
-  Z2 .= 0
-  @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))])
-  @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))])
-
-  # Step 4.3: Compute 
-  # Z₂ = I + Z₂
-  for i = 1:(2*p)
-    Z2[i, i] += 1
-  end
+    @views Z1[1:n, 1:p] .= Uk .* (-1)
+    @views Z1[1:n, (p+1):(2*p)] .= Vk
+    try
+      ma57_solve!(workspace.M, Z1, workspace._qn_work)
+    catch e
+      !(e isa HSL.Ma57Exception) && rethrow(e)
+    end
+    if any(isnan, Z1) || workspace.M.info.info[1] < 0
+      workspace.status = :failed
+      return
+    end
 
-  # Step 5: Solve
-  # (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )⁻¹[y₁]
-  # (       [A     -αI]     )  [y₁]
-  # using Julia LinearALgebra's lu!
-  F = lu!(Z2[1:(2*p), 1:(2*p)], check = false) # FIXME ?
-  @views ldiv!(y2[1:(2*p)], F, y1[1:(2*p)])
-  if any(isnan, y2)
-    workspace.status = :failed
-    return
-  end
+    # Step 4.2: Compute 
+    # Z₂ = FᵀZ₁ = UᵀZ₁[1:n]
+    # Z₂ = FᵀZ₁ = VᵀZ₁[1:n]
+    Z2 .= 0
+    @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))])
+    @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))])
+
+    # Step 4.3: Compute 
+    # Z₂ = I + Z₂
+    for i = 1:(2*p)
+      Z2[i, i] += 1
+    end
+
+    # Step 5: Solve
+    # (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )⁻¹[y₁]
+    # (       [A     -αI]     )  [y₁]
+    # using LAPACK
+    @views info_f =
+      getrf!(BlasInt(2p), BlasInt(2p), Z2, stride(Z2, 2), workspace._ipiv, workspace._info)
+    if info_f != 0
+      workspace.status = :failed
+      return
+    end
 
-  # Step 6: Compute
-  # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂]
-  # x₂ = E[y₂] = [ 0 0][y₂] = [0]
-  @views mul!(x2[1:n], Vk, y2[(p+1):(2*p)])
-  @views mul!(x2[1:n], Uk, y2[1:p], -one(eltype(y2)), one(eltype(y2)))
+    @views info_s = getrs!(
+      'N',
+      BlasInt(2p),
+      BlasInt(1),
+      Z2[1:(2p), 1:(2p)],
+      stride(Z2, 2),
+      workspace._ipiv,
+      y1[1:(2p)],
+      BlasInt(2p),
+      workspace._info,
+    )
+    if info_s != 0
+      workspace.status = :failed
+      return
+    end
+    if any(isnan, @view y1[1:(2*p)])
+      workspace.status = :failed
+      return
+    end
 
-  # Step 7: Solve
-  # [x₃] = [σI+ξI  Aᵀ]⁻¹[x₂]
-  # [x₃] = [A     -αI]  [x₂]
-  try
-    ma57_solve!(workspace.M, x2, x3, workspace.dx, workspace.work, 10)
-  catch e
-    !(e isa HSL.Ma57Exception) && rethrow(e)
-  end
-  if any(isnan, x3) || workspace.M.info.info[1] < 0
-    workspace.status = :failed
-    return
-  end
+    # Step 6: Compute
+    # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂]
+    # x₂ = E[y₂] = [ 0 0][y₂] = [0]
+    @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)])
+    @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1)))
 
-  # Step 8:
-  # [B  Aᵀ]⁻¹[u] = x₁ - x₃ 
-  # [A -αI]  [u] = x₁ - x₃
-  workspace.x .= x1 .- x3
+    # Step 7: Solve
+    # [x₃] = [σI+ξI  Aᵀ]⁻¹[x₂]
+    # [x₃] = [A     -αI]  [x₂]
+    try
+      ma57_solve!(workspace.M, x2, x3, workspace.dx, workspace.work, 10)
+    catch e
+      !(e isa HSL.Ma57Exception) && rethrow(e)
+    end
+    if any(isnan, x3) || workspace.M.info.info[1] < 0
+      workspace.status = :failed
+      return
+    end
+
+    # Step 8:
+    # [B  Aᵀ]⁻¹[u] = x₁ - x₃ 
+    # [A -αI]  [u] = x₁ - x₃
+    workspace.x .= x1 .- x3
+  else
+    workspace.x .= x1
+  end
 end
 
 function get_solution!(
diff --git a/ext/LDLFactorizations/ldlt.jl b/ext/LDLFactorizations/ldlt.jl
index 3cd46e1d..8305073f 100644
--- a/ext/LDLFactorizations/ldlt.jl
+++ b/ext/LDLFactorizations/ldlt.jl
@@ -2,6 +2,7 @@ mutable struct PenaltyLDLTWorkspace{
   WP<:LDLFactorization,
   K2<:AbstractMatrix,
   V<:AbstractVector,
+  VI<:Union{Nothing,AbstractVector},
   T<:Real,
 } <: AbstractLDLTWorkspace
   M::WP
@@ -9,6 +10,8 @@ mutable struct PenaltyLDLTWorkspace{
   x::V
   dx::V
   r::V
+  _ipiv::VI # For CompactBFGS LU factorization
+  _info::Base.RefValue{BlasInt} # For CompactBFGS LU factorization
   σ::T
   n::Int
   m::Int
@@ -39,6 +42,8 @@ function construct_ldlt_workspace(
     similar(u1),
     similar(u1),
     similar(u1),
+    nothing,
+    Ref{BlasInt}(0),
     zero(T),
     n,
     m,
@@ -60,6 +65,8 @@ function construct_ldlt_workspace(
     similar(u1),
     similar(u1),
     similar(u1),
+    Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem),
+    Ref{BlasInt}(0),
     zero(T),
     n,
     m,
@@ -268,12 +275,12 @@ function solve_system!(
   B = workspace.H.B
   n, m = workspace.n, workspace.m
   p = min(B._insert - 1, B._mem)
-  x1, x2, x3, y1, y2 = H.x1, H.x2, H.x3, H.y1, H.y2
+  x1, x2, x3, y1 = H.x1, H.x2, H.x3, H.y1
   Z1, Z2 = H.Z1, H.Z2
   Uk = @view B.Uk[:, 1:p]
   Vk = @view B.Vk[:, 1:p]
 
-  # Step 0: Write (#TODO: we can use easily use QRMumps instead of LDLFactorization here...)
+  # Step 0: Write
   # [B  Aᵀ] = [σI+ξI  Aᵀ] + [-U V]([U V])ᵀ
   # [A -αI] = [A     -αI] + [ 0 0]([0 0])
   # Hence,
@@ -292,7 +299,7 @@ function solve_system!(
     return
   end
 
-  # Step 2: Compute # TODO: allow for iterative refinement
+  # Step 2: Compute
   # [x₁] = [σI+ξI  Aᵀ]⁻¹[u]
   # [x₁] = [A     -αI]  [u]
   ldiv!(x1, workspace.M, u)
@@ -301,71 +308,96 @@ function solve_system!(
     return
   end
 
-  # Step 3: Compute
-  # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)]
-  # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)]
-  @views mul!(y1[1:p], Uk', x1[1:n])
-  @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n])
+  if p > 0
 
+    # Step 3: Compute
+    # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)]
+    # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)]
+    @views mul!(y1[1:p], Uk', x1[1:n])
+    @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n])
 
-  # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )
-  #                                   (       [A     -αI]     )
-  # Step 4.1: Compute 
-  # Z₁ = [σI+ξI  Aᵀ]⁻¹ E = [σI+ξI  Aᵀ]⁻¹[-U V]
-  # Z₁ = [A     -αI]   E = [A     -αI]  [ 0 0]
-  Z1 .= 0
 
-  @views Z1[1:n, 1:p] .= Uk .* (-1)
-  @views Z1[1:n, (p+1):(2*p)] .= Vk
-  ldiv!(workspace.M, Z1)
-  if any(isnan, Z1)
-    workspace.status = :failed
-    return
-  end
+    # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )
+    #                                   (       [A     -αI]     )
+    # Step 4.1: Compute 
+    # Z₁ = [σI+ξI  Aᵀ]⁻¹ E = [σI+ξI  Aᵀ]⁻¹[-U V]
+    # Z₁ = [A     -αI]   E = [A     -αI]  [ 0 0]
+    Z1 .= 0
 
-  # Step 4.2: Compute 
-  # Z₂ = FᵀZ₁ = UᵀZ₁[1:n]
-  # Z₂ = FᵀZ₁ = VᵀZ₁[1:n]
-  Z2 .= 0
-  @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))])
-  @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))])
-
-  # Step 4.3: Compute 
-  # Z₂ = I + Z₂
-  for i = 1:(2*p)
-    Z2[i, i] += 1
-  end
+    @views Z1[1:n, 1:p] .= Uk .* (-1)
+    @views Z1[1:n, (p+1):(2*p)] .= Vk
+    ldiv!(workspace.M, Z1)
+    if any(isnan, Z1)
+      workspace.status = :failed
+      return
+    end
 
-  # Step 5: Solve
-  # (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )⁻¹[y₁]
-  # (       [A     -αI]     )  [y₁]
-  # using Julia LinearALgebra's lu!
-  F = lu!(Z2[1:(2*p), 1:(2*p)], check = false) # FIXME ?
-  @views ldiv!(y2[1:(2*p)], F, y1[1:(2*p)])
-  if any(isnan, y2)
-    workspace.status = :failed
-    return
-  end
+    # Step 4.2: Compute 
+    # Z₂ = FᵀZ₁ = UᵀZ₁[1:n]
+    # Z₂ = FᵀZ₁ = VᵀZ₁[1:n]
+    Z2 .= 0
+    @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))])
+    @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))])
+
+    # Step 4.3: Compute 
+    # Z₂ = I + Z₂
+    for i = 1:(2*p)
+      Z2[i, i] += 1
+    end
 
-  # Step 6: Compute
-  # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂]
-  # x₂ = E[y₂] = [ 0 0][y₂] = [0]
-  @views mul!(x2[1:n], Vk, y2[(p+1):(2*p)])
-  @views mul!(x2[1:n], Uk, y2[1:p], -one(eltype(y2)), one(eltype(y2)))
-
-  # Step 7: Solve
-  # [x₃] = [σI+ξI  Aᵀ]⁻¹[x₂]
-  # [x₃] = [A     -αI]  [x₂]
-  ldiv!(x3, workspace.M, x2)
-  if any(isnan, x3)
-    workspace.status = :failed
-    return
-  end
+    # Step 5: Solve
+    # (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )⁻¹[y₁]
+    # (       [A     -αI]     )  [y₁]
+    # using LAPACK
+    @views info_f =
+      getrf!(BlasInt(2p), BlasInt(2p), Z2, stride(Z2, 2), workspace._ipiv, workspace._info)
+    if info_f != 0
+      workspace.status = :failed
+      return
+    end
 
-  # Step 8:
-  # [B  Aᵀ]⁻¹[u] = x₁ - x₃ 
-  # [A -αI]  [u] = x₁ - x₃
-  workspace.x .= x1 .- x3
+    @views info_s = getrs!(
+      'N',
+      BlasInt(2p),
+      BlasInt(1),
+      Z2[1:(2p), 1:(2p)],
+      stride(Z2, 2),
+      workspace._ipiv,
+      y1[1:(2p)],
+      BlasInt(2p),
+      workspace._info,
+    )
+    if info_s != 0
+      workspace.status = :failed
+      return
+    end
+    if any(isnan, @view y1[1:(2*p)])
+      workspace.status = :failed
+      return
+    end
+
+    # Step 6: Compute
+    # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂]
+    # x₂ = E[y₂] = [ 0 0][y₂] = [0]
+    @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)])
+    @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1)))
+
+    # Step 7: Solve
+    # [x₃] = [σI+ξI  Aᵀ]⁻¹[x₂]
+    # [x₃] = [A     -αI]  [x₂]
+    ldiv!(x3, workspace.M, x2)
+    if any(isnan, x3)
+      workspace.status = :failed
+      return
+    end
+
+    # Step 8:
+    # [B  Aᵀ]⁻¹[u] = x₁ - x₃ 
+    # [A -αI]  [u] = x₁ - x₃
+    workspace.x .= x1 .- x3
+  else
+    workspace.x .= x1
+  end
 end
 
 function get_solution!(x::V, workspace::PenaltyLDLTWorkspace) where {V<:AbstractVector}
diff --git a/ext/PeneloptHSLExt.jl b/ext/PeneloptHSLExt.jl
index b6c36327..5c6f5dfb 100644
--- a/ext/PeneloptHSLExt.jl
+++ b/ext/PeneloptHSLExt.jl
@@ -5,9 +5,14 @@ using Penelopt
 
 using LinearAlgebra, SparseMatricesCOO
 
+# Import BLAS functions
+import LinearAlgebra.BLAS: @blasfunc
+import LinearAlgebra: BlasInt, libblastrampoline
+
 import Penelopt: AbstractHSLWorkspace
 import Penelopt: construct_ma57_workspace, solve_system!, update_workspace!
 import Penelopt: get_inertia, get_solution!, get_status
+import Penelopt: getrf!, getrs!
 import Penelopt: set_dual_inertia!, set_primal_inertia!
 
 function __init__()
diff --git a/ext/PeneloptLDLFactorizationsExt.jl b/ext/PeneloptLDLFactorizationsExt.jl
index c30ab3ad..d8990bb8 100644
--- a/ext/PeneloptLDLFactorizationsExt.jl
+++ b/ext/PeneloptLDLFactorizationsExt.jl
@@ -5,9 +5,14 @@ using Penelopt
 
 using LinearAlgebra, SparseArrays, SparseMatricesCOO
 
+# Import BLAS functions
+import LinearAlgebra.BLAS: @blasfunc
+import LinearAlgebra: BlasInt, libblastrampoline
+
 import Penelopt: AbstractLDLTWorkspace
 import Penelopt: construct_ldlt_workspace, solve_system!, update_workspace!
 import Penelopt: get_inertia, get_solution!, get_status
+import Penelopt: getrf!, getrs!
 import Penelopt: set_dual_inertia!, set_primal_inertia!
 
 include("LDLFactorizations/ldlt.jl")
diff --git a/src/Penelopt.jl b/src/Penelopt.jl
index c5aac526..8c1c1d66 100644
--- a/src/Penelopt.jl
+++ b/src/Penelopt.jl
@@ -22,6 +22,10 @@ using NLPModels, NLPModelsModifiers
 using LinearOperators, QuadraticModels, SolverCore, SparseMatricesCOO
 using MPI, MUMPS
 
+# Import BLAS functions
+import LinearAlgebra.BLAS: @blasfunc
+import LinearAlgebra: BlasInt, libblastrampoline
+
 import SolverCore: get_status, reset!
 
 function __init__()
@@ -47,6 +51,7 @@ include("types/pre-processing/FixedVariable.jl")
 include("linear_algebra/K2.jl")
 include("linear_algebra/construct_workspace.jl")
 include("linear_algebra/mumps.jl")
+include("linear_algebra/lapack.jl")
 
 include("types/PenalizedProblem.jl")
 include("types/ShiftedPenalizedProblem.jl")
diff --git a/src/linear_algebra/K2.jl b/src/linear_algebra/K2.jl
index c7e84255..2381d285 100644
--- a/src/linear_algebra/K2.jl
+++ b/src/linear_algebra/K2.jl
@@ -15,7 +15,6 @@ mutable struct CompactBFGSK2{
   x2::V
   x3::V
   y1::V
-  y2::V
 end
 
 function K2(
@@ -128,7 +127,6 @@ function K2(
     zeros(T, n+m),
     zeros(T, n+m),
     zeros(T, 2*B._mem),
-    zeros(T, 2*B._mem),
   )
 end
 
diff --git a/src/linear_algebra/lapack.jl b/src/linear_algebra/lapack.jl
new file mode 100644
index 00000000..334643dd
--- /dev/null
+++ b/src/linear_algebra/lapack.jl
@@ -0,0 +1,66 @@
+for (getrf, getrs, T) in ((:sgetrf_, :sgetrs_, :Float32), (:dgetrf_, :dgetrs_, :Float64))
+  @eval begin
+
+    function getrf!(
+      m::BlasInt,
+      n::BlasInt,
+      a::AbstractMatrix{$T},
+      lda::BlasInt,
+      ipiv::AbstractVector{BlasInt},
+      info::Base.RefValue{BlasInt},
+    )
+      ccall(
+        (@blasfunc($getrf), libblastrampoline),
+        Cvoid,
+        (Ref{BlasInt}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}),
+        m,
+        n,
+        a,
+        lda,
+        ipiv,
+        info,
+      )
+      return info[]
+    end
+
+    function getrs!(
+      trans::Char,
+      n::BlasInt,
+      nrhs::BlasInt,
+      a::AbstractMatrix{$T},
+      lda::BlasInt,
+      ipiv::AbstractVector{BlasInt},
+      b::AbstractVecOrMat{$T},
+      ldb::BlasInt,
+      info::Base.RefValue{BlasInt},
+    )
+      ccall(
+        (@blasfunc($getrs), libblastrampoline),
+        Cvoid,
+        (
+          Ref{UInt8},
+          Ref{BlasInt},
+          Ref{BlasInt},
+          Ptr{$T},
+          Ref{BlasInt},
+          Ptr{BlasInt},
+          Ptr{$T},
+          Ref{BlasInt},
+          Ref{BlasInt},
+          Clong,
+        ),
+        trans,
+        n,
+        nrhs,
+        a,
+        lda,
+        ipiv,
+        b,
+        ldb,
+        info,
+        1,
+      )
+      return info[]
+    end
+  end
+end
diff --git a/src/linear_algebra/mumps.jl b/src/linear_algebra/mumps.jl
index c98e6f4a..4f0c3422 100644
--- a/src/linear_algebra/mumps.jl
+++ b/src/linear_algebra/mumps.jl
@@ -2,11 +2,14 @@ mutable struct PenaltyMUMPSWorkspace{
   WP<:Mumps,
   K2<:AbstractMatrix,
   V<:AbstractVector,
+  VI<:Union{Nothing,AbstractVector},
   T<:Real,
 } <: AbstractMUMPSWorkspace
   M::WP
   H::K2
   x::V
+  _ipiv::VI # For CompactBFGS LU factorization
+  _info::Base.RefValue{BlasInt}  # For CompactBFGS LU factorization
   σ::T
   n::Int
   m::Int
@@ -84,7 +87,19 @@ function construct_mumps_workspace(
   S.rhs = pointer(x)
   S._y_gc_haven = x
 
-  return PenaltyMUMPSWorkspace(S, H, x, zero(T), n, m, :uninitialized, false, 0)
+  return PenaltyMUMPSWorkspace(
+    S,
+    H,
+    x,
+    nothing,
+    Ref{BlasInt}(0),
+    zero(T),
+    n,
+    m,
+    :uninitialized,
+    false,
+    0,
+  )
 end
 
 function construct_mumps_workspace(
@@ -145,7 +160,19 @@ function construct_mumps_workspace(
   S.rhs = pointer(x)
   S._y_gc_haven = x
 
-  return PenaltyMUMPSWorkspace(S, H, x, zero(T), n, m, :uninitialized, false, 0)
+  return PenaltyMUMPSWorkspace(
+    S,
+    H,
+    x,
+    Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem),
+    Ref{BlasInt}(0),
+    zero(T),
+    n,
+    m,
+    :uninitialized,
+    false,
+    0,
+  )
 end
 
 function update_workspace!(
@@ -280,7 +307,7 @@ function solve_system!(
   mumps = workspace.M
   n, m = workspace.n, workspace.m
   p = min(B._insert - 1, B._mem)
-  x1, x2, x3, y1, y2 = H.x1, H.x2, H.x3, H.y1, H.y2
+  x1, x2, x3, y1 = H.x1, H.x2, H.x3, H.y1
   Z1, Z2 = H.Z1, H.Z2
 
   Uk = @view B.Uk[:, 1:p]
@@ -332,7 +359,7 @@ function solve_system!(
   # [x₁] = [A     -αI]  [u]
 
   x1 .= u
-  MUMPS.associate_rhs!(mumps, x1)
+  MUMPS.associate_rhs!(mumps, x1; unsafe = true)
   MUMPS.mumps_solve!(x1, mumps; rhs_changed = true)
 
   # MUMPS infog(1): a negative value is an error in the factorization.
@@ -343,83 +370,108 @@ function solve_system!(
 
   update_pivtol!(workspace)
 
-  # Step 3: Compute
-  # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)]
-  # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)]
-  @views mul!(y1[1:p], Uk', x1[1:n])
-  @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n])
+  if p > 0
 
+    # Step 3: Compute
+    # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)]
+    # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)]
+    @views mul!(y1[1:p], Uk', x1[1:n])
+    @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n])
 
-  # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )
-  #                                   (       [A     -αI]     )
-  # Step 4.1: Compute 
-  # Z₁ = [σI+ξI  Aᵀ]⁻¹ E = [σI+ξI  Aᵀ]⁻¹[-U V]
-  # Z₁ = [A     -αI]   E = [A     -αI]  [ 0 0]
-  Z1 .= 0
 
-  @views Z1[1:n, 1:p] .= Uk .* (-1)
-  @views Z1[1:n, (p+1):(2*p)] .= Vk
+    # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )
+    #                                   (       [A     -αI]     )
+    # Step 4.1: Compute 
+    # Z₁ = [σI+ξI  Aᵀ]⁻¹ E = [σI+ξI  Aᵀ]⁻¹[-U V]
+    # Z₁ = [A     -αI]   E = [A     -αI]  [ 0 0]
+    Z1 .= 0
 
-  MUMPS.associate_rhs!(mumps, Z1)
-  MUMPS.mumps_solve!(Z1, mumps; rhs_changed = true)
+    @views Z1[1:n, 1:p] .= Uk .* (-1)
+    @views Z1[1:n, (p+1):(2*p)] .= Vk
 
-  # MUMPS infog(1): a negative value is an error in the factorization.
-  if any(isnan, Z1) || mumps.infog[1] < 0
-    workspace.status = :failed
-    return
-  end
+    MUMPS.associate_rhs!(mumps, Z1; unsafe = true)
+    MUMPS.mumps_solve!(Z1, mumps; rhs_changed = true)
 
-  update_pivtol!(workspace)
+    # MUMPS infog(1): a negative value is an error in the factorization.
+    if any(isnan, Z1) || mumps.infog[1] < 0
+      workspace.status = :failed
+      return
+    end
 
-  # Step 4.2: Compute 
-  # Z₂ = FᵀZ₁ = UᵀZ₁[1:n]
-  # Z₂ = FᵀZ₁ = VᵀZ₁[1:n]
-  Z2 .= 0
-  @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))])
-  @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))])
-
-  # Step 4.3: Compute 
-  # Z₂ = I + Z₂
-  for i = 1:(2*p)
-    Z2[i, i] += 1
-  end
+    update_pivtol!(workspace)
 
-  # Step 5: Solve
-  # (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )⁻¹[y₁]
-  # (       [A     -αI]     )  [y₁]
-  # using Julia LinearALgebra's lu!
-  F = lu!(Z2[1:(2*p), 1:(2*p)], check = false) # FIXME ?
-  @views ldiv!(y2[1:(2*p)], F, y1[1:(2*p)])
-  if any(isnan, y2)
-    workspace.status = :failed
-    return
-  end
+    # Step 4.2: Compute 
+    # Z₂ = FᵀZ₁ = UᵀZ₁[1:n]
+    # Z₂ = FᵀZ₁ = VᵀZ₁[1:n]
+    Z2 .= 0
+    @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))])
+    @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))])
 
-  # Step 6: Compute
-  # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂]
-  # x₂ = E[y₂] = [ 0 0][y₂] = [0]
-  @views mul!(x2[1:n], Vk, y2[(p+1):(2*p)])
-  @views mul!(x2[1:n], Uk, y2[1:p], -one(eltype(y2)), one(eltype(y2)))
+    # Step 4.3: Compute 
+    # Z₂ = I + Z₂
+    for i = 1:(2*p)
+      Z2[i, i] += 1
+    end
 
-  # Step 7: Solve
-  # [x₃] = [σI+ξI  Aᵀ]⁻¹[x₂]
-  # [x₃] = [A     -αI]  [x₂]
-  x3 .= x2
-  MUMPS.associate_rhs!(mumps, x3)
-  MUMPS.mumps_solve!(x3, mumps; rhs_changed = true)
+    # Step 5: Solve
+    # (I + Fᵀ [σI+ξI  Aᵀ]⁻¹ E )⁻¹[y₁]
+    # (       [A     -αI]     )  [y₁]
+    # using LAPACK
+    @views info_f =
+      getrf!(BlasInt(2p), BlasInt(2p), Z2, stride(Z2, 2), workspace._ipiv, workspace._info)
+    if info_f != 0
+      workspace.status = :failed
+      return
+    end
 
-  # MUMPS infog(1): a negative value is an error in the factorization.
-  if any(isnan, x3) || mumps.infog[1] < 0
-    workspace.status = :failed
-    return
-  end
+    @views info_s = getrs!(
+      'N',
+      BlasInt(2p),
+      BlasInt(1),
+      Z2[1:(2p), 1:(2p)],
+      stride(Z2, 2),
+      workspace._ipiv,
+      y1[1:(2p)],
+      BlasInt(2p),
+      workspace._info,
+    )
+    if info_s != 0
+      workspace.status = :failed
+      return
+    end
+    if any(isnan, @view y1[1:(2*p)])
+      workspace.status = :failed
+      return
+    end
 
-  update_pivtol!(workspace)
+    # Step 6: Compute
+    # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂]
+    # x₂ = E[y₂] = [ 0 0][y₂] = [0]
+    @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)])
+    @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1)))
+
+    # Step 7: Solve
+    # [x₃] = [σI+ξI  Aᵀ]⁻¹[x₂]
+    # [x₃] = [A     -αI]  [x₂]
+    x3 .= x2
+    MUMPS.associate_rhs!(mumps, x3; unsafe = true)
+    MUMPS.mumps_solve!(x3, mumps; rhs_changed = true)
 
-  # Step 8:
-  # [B  Aᵀ]⁻¹[u] = x₁ - x₃ 
-  # [A -αI]  [u] = x₁ - x₃
-  workspace.x .= x1 .- x3
+    # MUMPS infog(1): a negative value is an error in the factorization.
+    if any(isnan, x3) || mumps.infog[1] < 0
+      workspace.status = :failed
+      return
+    end
+
+    update_pivtol!(workspace)
+
+    # Step 8:
+    # [B  Aᵀ]⁻¹[u] = x₁ - x₃ 
+    # [A -αI]  [u] = x₁ - x₃
+    workspace.x .= x1 .- x3
+  else
+    workspace.x .= x1
+  end
 end
 
 function get_solution!(x::V, workspace::PenaltyMUMPSWorkspace) where {V<:AbstractVector}
diff --git a/src/types/PenalizedProblem.jl b/src/types/PenalizedProblem.jl
index bbcef241..da75f830 100644
--- a/src/types/PenalizedProblem.jl
+++ b/src/types/PenalizedProblem.jl
@@ -37,14 +37,12 @@ function L2PenalizedProblem(nlp::AbstractNLPModel{T,S}) where {T,S}
   A = SparseMatrixCOO(nlp.meta.ncon, nlp.meta.nvar, rows, cols, vals)
   b = similar(x0, eltype(x0), nlp.meta.ncon)
 
-  store_previous_jacobian = isa(nlp, QuasiNewtonModel) ? true : false
   penalty = CompositeNormL2(
     one(T),
     (c, x) -> cons!(nlp, x, c),
     (j, x) -> jac_coord!(nlp, x, j.vals),
     A,
     b,
-    store_previous_jacobian = store_previous_jacobian,
   )
   return L2PenalizedProblem(nlp, penalty, nlp.meta)
 end
diff --git a/src/types/ShiftedPenalizedProblem.jl b/src/types/ShiftedPenalizedProblem.jl
index d3dc0779..7e17c99d 100644
--- a/src/types/ShiftedPenalizedProblem.jl
+++ b/src/types/ShiftedPenalizedProblem.jl
@@ -31,7 +31,6 @@ mutable struct ShiftedL2PenalizedProblem{
   h::H
   parent::P
   meta::meta
-  _qn_∇f_prev::SN
   _qn_y::SN
   _qn_x_prev::SN
   _is_first_shift::Bool
@@ -56,7 +55,6 @@ function ShiftedL2PenalizedProblem(
     ψ,
     penalty_nlp,
     penalty_nlp.meta,
-    zero(∇f),
     similar(∇f),
     zero(∇f),
     true,
@@ -91,7 +89,6 @@ function ShiftedL2PenalizedProblem(
     penalty_nlp.meta,
     nothing,
     nothing,
-    nothing,
     true,
   )
 end
@@ -163,23 +160,27 @@ function shift!(
 }
   nlp, h = shifted_penalty_nlp.parent.model, shifted_penalty_nlp.parent.h
   φ, ψ = shifted_penalty_nlp.model, shifted_penalty_nlp.h
-  qn_y, qn_g_prev, qn_x_prev = shifted_penalty_nlp._qn_y,
-  shifted_penalty_nlp._qn_∇f_prev,
-  shifted_penalty_nlp._qn_x_prev
+  qn_y, qn_x_prev = shifted_penalty_nlp._qn_y, shifted_penalty_nlp._qn_x_prev
   is_first_shift = shifted_penalty_nlp._is_first_shift
 
   qn_s = qn_x_prev
   g, B = φ.data.c, φ.data.H
 
+  if !is_first_shift
+    qn_y .= g
+    if !isnothing(y)
+      mul!(qn_y, ψ.A', y, -one(T), -one(T)) # y = - g_prev - J(x)_prev^T λ
+    end
+  end
+
   isnothing(∇f) ? grad!(nlp, x, g) : (g .= ∇f)
   shift!(ψ, x, J = J, c = c)
 
   # Update the approximation.
   if !is_first_shift
-    @. qn_y = g - qn_g_prev
+    @. qn_y .+= g
     if !isnothing(y)
       mul!(qn_y, ψ.A', y, one(T), one(T)) # y = y + J(x)^T λ 
-      mul!(qn_y, ψ.A_prev', y, -one(T), one(T)) # y = y - J(x)_prev^T λ
     end
 
     qn_s .= x .- qn_x_prev
@@ -189,11 +190,7 @@ function shift!(
     shifted_penalty_nlp._is_first_shift = false
   end
 
-  # Copy the gradient and Jacobian.
-  qn_g_prev .= g
-  ψ.A_prev.vals .= ψ.A.vals
   qn_x_prev .= x
-
 end
 
 function shift!(
@@ -254,9 +251,7 @@ function reset!(
   nlp, h = shifted_penalty_nlp.parent.model, shifted_penalty_nlp.parent.h
   φ, ψ = shifted_penalty_nlp.model, shifted_penalty_nlp.h
   x_prev = shifted_penalty_nlp._qn_x_prev .= 0
-  g_prev = shifted_penalty_nlp._qn_∇f_prev .= 0
 
-  ψ.A_prev.vals .= ψ.A.vals
   LinearOperators.reset!(φ.data.H)
   shifted_penalty_nlp._is_first_shift = true
 end
diff --git a/src/types/norm/CompositeNormL2.jl b/src/types/norm/CompositeNormL2.jl
index 12a45949..4c96dea4 100644
--- a/src/types/norm/CompositeNormL2.jl
+++ b/src/types/norm/CompositeNormL2.jl
@@ -1,7 +1,7 @@
 abstract type AbstractCompositeNorm end
 
 @doc raw"""
-    CompositeNormL2(λ, c!, J!, A, b; store_previous_jacobian::Bool = false)
+    CompositeNormL2(λ, c!, J!, A, b)
 
 Returns function `c` composed with the `ℓ₂` norm:
 ```math
@@ -21,8 +21,6 @@ such that `J` is the Jacobian of `c`. It is expected that `m ≤ n`.
     c!(b <: AbstractVector{Real}, xk <: AbstractVector{Real})
     J!(A <: AbstractSparseMatrixCOO{Real, Integer}, xk <: AbstractVector{Real})
 ```
-Moreover, if you want shifted instances of the operator to store the previous Jacobian on each shift, you can specify `store_previous_jacobian = true`.
-This is particularly useful for quasi-Newton updates in the context of constrained optimization.
 """
 mutable struct CompositeNormL2{
   T<:Real,
@@ -36,7 +34,6 @@ mutable struct CompositeNormL2{
   J!::F1
   A::M
   b::V
-  store_previous_jacobian::Bool
 
   function CompositeNormL2(
     λ::T,
@@ -44,20 +41,12 @@ mutable struct CompositeNormL2{
     J!::Function,
     A::AbstractMatrix{T},
     b::AbstractVector{T};
-    store_previous_jacobian::Bool = false,
   ) where {T<:Real}
     λ > 0 || error("CompositeNormL2: λ should be positive")
     length(b) == size(A, 1) || error(
       "Composite Norm L2: Wrong input dimensions, the length of c(x) should be the same as the number of rows of J(x)",
     )
-    new{T,typeof(c!),typeof(J!),typeof(A),typeof(b)}(
-      NormL2(λ),
-      c!,
-      J!,
-      A,
-      b,
-      store_previous_jacobian,
-    )
+    new{T,typeof(c!),typeof(J!),typeof(A),typeof(b)}(NormL2(λ), c!, J!, A, b)
   end
 end
 
diff --git a/src/types/norm/ShiftedCompositeNormL2.jl b/src/types/norm/ShiftedCompositeNormL2.jl
index 4b3126ad..2f6afa32 100644
--- a/src/types/norm/ShiftedCompositeNormL2.jl
+++ b/src/types/norm/ShiftedCompositeNormL2.jl
@@ -1,7 +1,7 @@
 abstract type AbstractShiftedCompositeNorm end
 
 @doc raw"""
-    ShiftedCompositeNormL2(h, c!, J!, A, b; store_previous_jacobian::Bool = false)
+    ShiftedCompositeNormL2(h, c!, J!, A, b)
 
 Returns the shift of a function `c` composed with the `ℓ₂` norm (see CompositeNormL2.jl).
 Here, `c` is linearized i.e, `c(x + s) ≈ c(x) + J(x)s`. 
@@ -22,23 +22,18 @@ such that `J` is the Jacobian of `c`. It is expected that `m ≤ n`.
 c!(b <: AbstractVector{Real}, xk <: AbstractVector{Real})
 J!(A <: AbstractSparseMatrixCOO{Real, Integer}, xk <: AbstractVector{Real})
 ```
-Moreover, if you want shifted instances of the operator to store the previous Jacobian on each shift, you can specify `store_previous_jacobian = true`.
-In this case, each time a shift is performed, the previous Jacobian is stored in the `A_prev` field.
-This is particularly useful for quasi-Newton updates in the context of constrained optimization.
 """
 mutable struct ShiftedCompositeNormL2{
   T<:Real,
   F0<:Function,
   F1<:Function,
   M<:AbstractMatrix{T},
-  N<:Union{Nothing,M},
   V<:AbstractVector{T},
 } <: AbstractShiftedCompositeNorm
   h::NormL2{T}
   c!::F0
   J!::F1
   A::M
-  A_prev::N # (Optional) can be used to store the previous Jacobian, useful for quasi-Newton approximations
   b::V
   g::V
   function ShiftedCompositeNormL2(
@@ -47,7 +42,6 @@ mutable struct ShiftedCompositeNormL2{
     J!::Function,
     A::AbstractMatrix{T},
     b::AbstractVector{T};
-    store_previous_jacobian::Bool = false,
   ) where {T<:Real}
     if length(b) != size(A, 1)
       error(
@@ -55,18 +49,9 @@ mutable struct ShiftedCompositeNormL2{
       )
     end
 
-    A_prev = store_previous_jacobian ? copy(A) : nothing
     g = similar(b)
 
-    new{T,typeof(c!),typeof(J!),typeof(A),typeof(A_prev),typeof(b)}(
-      NormL2(λ),
-      c!,
-      J!,
-      A,
-      A_prev,
-      b,
-      g,
-    )
+    new{T,typeof(c!),typeof(J!),typeof(A),typeof(b)}(NormL2(λ), c!, J!, A, b, g)
   end
 end
 
@@ -77,14 +62,7 @@ shifted(
   begin
     A = ψ.A
     b = similar(ψ.b)
-    ShiftedCompositeNormL2(
-      ψ.h.lambda,
-      ψ.c!,
-      ψ.J!,
-      A,
-      b,
-      store_previous_jacobian = ψ.store_previous_jacobian,
-    )
+    ShiftedCompositeNormL2(ψ.h.lambda, ψ.c!, ψ.J!, A, b)
   end
 
 fun_name(ψ::ShiftedCompositeNormL2) = "shifted `ℓ₂` norm"
@@ -103,7 +81,6 @@ function shift!(
   J = nothing,
   c = nothing,
 ) where {R<:Real}
-  !isnothing(ψ.A_prev) && (ψ.A_prev.vals .= ψ.A.vals) # Update previous Jacobian if necessary
   isnothing(c) ? ψ.c!(ψ.b, shift) : (ψ.b .= c)
   isnothing(J) ? ψ.J!(ψ.A, shift) : (ψ.A.vals .= J.vals)
   return ψ
diff --git a/test/instances/instance-generator.jl b/test/instances/instance-generator.jl
index 3c1d48a0..4ce88a11 100644
--- a/test/instances/instance-generator.jl
+++ b/test/instances/instance-generator.jl
@@ -77,15 +77,6 @@ function generate_instance(
     end
   end
   h = ShiftedCompositeNormL2(tau, c!, J!, SparseMatrixCOO(J), b)
-  return ShiftedL2PenalizedProblem(
-    model,
-    h,
-    nothing,
-    model.meta,
-    nothing,
-    nothing,
-    nothing,
-    true,
-  ),
+  return ShiftedL2PenalizedProblem(model, h, nothing, model.meta, nothing, nothing, true),
   Dict(:u => u, :y => y, :tau => tau)
 end
diff --git a/test/instances/instance-reader.jl b/test/instances/instance-reader.jl
index b52c5d83..11f0eecd 100644
--- a/test/instances/instance-reader.jl
+++ b/test/instances/instance-reader.jl
@@ -45,14 +45,5 @@ function read_instance(file::String; type = Float64, Hessian_modifier = H -> H)
     end
   end
   h = ShiftedCompositeNormL2(data["tau"], c!, J!, SparseMatrixCOO(data["J"]), data["c"])
-  return ShiftedL2PenalizedProblem(
-    model,
-    h,
-    nothing,
-    model.meta,
-    nothing,
-    nothing,
-    nothing,
-    true,
-  )
+  return ShiftedL2PenalizedProblem(model, h, nothing, model.meta, nothing, nothing, true)
 end
diff --git a/test/test-cutest.jl b/test/test-cutest.jl
index d4bfcf22..63ba7cc1 100644
--- a/test/test-cutest.jl
+++ b/test/test-cutest.jl
@@ -102,7 +102,13 @@ function test_problem(
 
     solver = L2PenaltySolver(LBFGS_model, linear_solver = linear_solver)
     stats_optimized = PeneloptExecutionStats(LBFGS_model)
-    solve!(solver, LBFGS_model, stats_optimized, atol = 1e-3, rtol = 1e-3, τ0 = 1.0)
+    if linear_solver != "mumps"
+      @test @wrappedallocs(
+        solve!(solver, LBFGS_model, stats_optimized, atol = 1e-3, rtol = 1e-3, τ0 = 1.0)
+      ) == 0
+    else
+      solve!(solver, LBFGS_model, stats_optimized, atol = 1e-3, rtol = 1e-3, τ0 = 1.0)
+    end
 
     stats_optimized.solution =
       recover_full_solution(LBFGS_model, stats_optimized.solution)