diff --git a/examples/envelope/envelope.jl b/examples/envelope/envelope.jl
index d83ab0343..d899bd2f0 100644
--- a/examples/envelope/envelope.jl
+++ b/examples/envelope/envelope.jl
@@ -67,7 +67,7 @@ function run_envelope()
 
     Hypatia.check_data(c, A, b, G, h, cone)
     (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true)
-    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, Q2, RiQ1)
+    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1)
 
     opt = Hypatia.Optimizer(maxiter=100, verbose=false)
     Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L)
diff --git a/examples/lp/lp.jl b/examples/lp/lp.jl
index 3bb8261eb..47c8c6bc9 100644
--- a/examples/lp/lp.jl
+++ b/examples/lp/lp.jl
@@ -60,7 +60,7 @@ function run_lp()
 
     Hypatia.check_data(c, A, b, G, h, cone)
     (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true)
-    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, Q2, RiQ1)
+    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1)
 
     opt = Hypatia.Optimizer(maxiter=100, verbose=false)
     Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L)
diff --git a/examples/namedpoly/namedpoly.jl b/examples/namedpoly/namedpoly.jl
index 050c3256b..bf368fa32 100644
--- a/examples/namedpoly/namedpoly.jl
+++ b/examples/namedpoly/namedpoly.jl
@@ -98,7 +98,7 @@ function run_namedpoly()
 
     Hypatia.check_data(c, A, b, G, h, cone)
     (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true)
-    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, Q2, RiQ1)
+    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1)
 
     opt = Hypatia.Optimizer(maxiter=100, verbose=false)
     Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L)
diff --git a/src/cone.jl b/src/cone.jl
index ef27536c0..bed5ba6af 100644
--- a/src/cone.jl
+++ b/src/cone.jl
@@ -46,11 +46,6 @@ function calcg!(g::Vector{Float64}, cone::Cone)
     return g
 end
 
-
-
-# TODO rewrite to return block-diagonal matrix type. maybe use
-# https://github.com/JuliaMatrices/BlockBandedMatrices.jl
-
 function calcHarr!(prod::AbstractMatrix{Float64}, arr::AbstractMatrix{Float64}, cone::Cone)
     for k in eachindex(cone.prms)
         calcHarr_prm!(view(prod, cone.idxs[k], :), view(arr, cone.idxs[k], :), cone.prms[k])
@@ -79,11 +74,6 @@ function calcHiarr!(prod::AbstractVector{Float64}, arr::AbstractVector{Float64},
     return prod
 end
 
-
-
-
-
-
 # utilities for converting between smat and svec forms (lower triangle) for symmetric matrices
 # TODO only need to do lower triangle if use symmetric matrix types
 const rt2 = sqrt(2)
diff --git a/src/linsyssolvers/naive.jl b/src/linsyssolvers/naive.jl
index 1ad529c46..7d9d751d5 100644
--- a/src/linsyssolvers/naive.jl
+++ b/src/linsyssolvers/naive.jl
@@ -3,6 +3,7 @@ naive method that simply performs one high-dimensional linear system solve
 TODO currently only does dense operations, needs to work for sparse
 =#
 mutable struct NaiveCache <: LinSysCache
+    cone
     c
     A
     b
@@ -26,10 +27,12 @@ mutable struct NaiveCache <: LinSysCache
         b::Vector{Float64},
         G::AbstractMatrix{Float64},
         h::Vector{Float64},
+        cone::Cone,
         )
 
         (n, p, q) = (length(c), length(b), length(h))
         L = new()
+        L.cone = cone
         L.c = c
         L.A = A
         L.b = b
@@ -69,8 +72,9 @@ function solvelinsys3!(
     rhs_tx::Vector{Float64},
     rhs_ty::Vector{Float64},
     rhs_tz::Vector{Float64},
-    H::AbstractMatrix{Float64},
-    L::NaiveCache,
+    mu::Float64,
+    L::NaiveCache;
+    identityH::Bool = false,
     )
 
     rhs = L.rhs3
@@ -79,7 +83,13 @@ function solvelinsys3!(
     @. rhs[L.tzk:L.tkk-1] = -rhs_tz
 
     @. L.LHS3copy = L.LHS3
-    @. L.LHS3copy[L.tzk:L.tkk-1, L.tzk:L.tkk-1] = -H
+    if !identityH
+        for k in eachindex(L.cone.prms)
+            idxs = L.tzk - 1 .+ L.cone.idxs[k]
+            dim = dimension(L.cone.prms[k])
+            calcHarr_prm!(view(L.LHS3copy, idxs, idxs), Matrix(-mu*I, dim, dim), L.cone.prms[k])
+        end
+    end
 
     F = bunchkaufman!(Symmetric(L.LHS3copy))
     ldiv!(F, rhs)
@@ -103,7 +113,6 @@ function solvelinsys6!(
     rhs_tau::Float64,
     mu::Float64,
     tau::Float64,
-    H::AbstractMatrix{Float64},
     L::NaiveCache,
     )
 
@@ -116,7 +125,11 @@ function solvelinsys6!(
     rhs[end] = rhs_tau
 
     @. L.LHS6copy = L.LHS6
-    L.LHS6copy[L.tzk:L.tkk-1, L.tsk:L.ttk-1] = H
+    for k in eachindex(L.cone.prms)
+        dim = dimension(L.cone.prms[k])
+        calcHarr_prm!(view(L.LHS6copy, L.tzk - 1 .+ L.cone.idxs[k], L.tsk - 1 .+ L.cone.idxs[k]), Matrix(mu*I, dim, dim), L.cone.prms[k])
+    end
+
     L.LHS6copy[L.tkk, end] = mu/tau/tau
 
     F = qr!(L.LHS6copy)
diff --git a/src/linsyssolvers/qrsymm.jl b/src/linsyssolvers/qrsymm.jl
index dde3098a4..b82b0f631 100644
--- a/src/linsyssolvers/qrsymm.jl
+++ b/src/linsyssolvers/qrsymm.jl
@@ -11,6 +11,7 @@ where M = -I (for initial iterate only) or M = -Hi (Hi is Hessian inverse, pre-s
 mutable struct QRSymmCache <: LinSysCache
     # TODO can remove some of the intermediary prealloced arrays after github.com/JuliaLang/julia/issues/23919 is resolved
     useiterative
+    cone
     c
     b
     G
@@ -42,6 +43,7 @@ mutable struct QRSymmCache <: LinSysCache
         b::Vector{Float64},
         G::AbstractMatrix{Float64},
         h::Vector{Float64},
+        cone::Cone,
         Q2::AbstractMatrix{Float64},
         RiQ1::AbstractMatrix{Float64};
         useiterative::Bool = false,
@@ -51,6 +53,7 @@ mutable struct QRSymmCache <: LinSysCache
         (n, p, q) = (length(c), length(b), length(h))
         nmp = n - p
         L.useiterative = useiterative
+        L.cone = cone
         L.c = c
         L.b = b
         L.G = G
@@ -87,7 +90,8 @@ QRSymmCache(
     A::AbstractMatrix{Float64},
     b::Vector{Float64},
     G::AbstractMatrix{Float64},
-    h::Vector{Float64};
+    h::Vector{Float64},
+    cone::Cone;
     useiterative::Bool = false,
     ) = error("to use a QRSymmCache for linear system solves, the data must be preprocessed and Q2 and RiQ1 must be passed into the QRSymmCache constructor")
 
@@ -96,11 +100,12 @@ function solvelinsys3!(
     rhs_tx::Vector{Float64},
     rhs_ty::Vector{Float64},
     rhs_tz::Vector{Float64},
-    H::AbstractMatrix{Float64},
-    L::QRSymmCache,
+    mu::Float64,
+    L::QRSymmCache;
+    identityH::Bool = false,
     )
 
-    F = helplhs!(H, L)
+    F = helplhs!(mu, L, identityH=identityH)
     @. rhs_ty *= -1.0
     @. rhs_tz *= -1.0
     helplinsys!(rhs_tx, rhs_ty, rhs_tz, F, L)
@@ -118,26 +123,26 @@ function solvelinsys6!(
     rhs_tau::Float64,
     mu::Float64,
     tau::Float64,
-    H::AbstractMatrix{Float64},
     L::QRSymmCache,
     )
 
     # solve two symmetric systems and combine the solutions
-    F = helplhs!(H, L)
+    F = helplhs!(mu, L)
 
     # (x2, y2, z2) = (rhs_tx, -rhs_ty, -H*rhs_ts - rhs_tz)
     @. rhs_ty *= -1.0
     @. rhs_tz *= -1.0
     if !iszero(rhs_ts)
-        mul!(L.z1, H, rhs_ts)
-        @. rhs_tz -= L.z1
+        calcHarr!(L.z1, rhs_ts, L.cone)
+        @. rhs_tz -= mu*L.z1
     end
     helplinsys!(rhs_tx, rhs_ty, rhs_tz, F, L)
 
     # (x1, y1, z1) = (-c, b, H*h)
     @. L.x1 = -L.c
     @. L.y1 = L.b
-    mul!(L.z1, H, L.h)
+    calcHarr!(L.z1, L.h, L.cone)
+    @. L.z1 *= mu
     helplinsys!(L.x1, L.y1, L.z1, F, L)
 
     # combine
@@ -196,12 +201,18 @@ end
 
 # calculate or factorize LHS of symmetric positive definite linear system
 function helplhs!(
-    H::AbstractMatrix{Float64},
-    L::QRSymmCache,
+    mu::Float64,
+    L::QRSymmCache;
+    identityH::Bool = false,
     )
 
     # Q2' * G' * H * G * Q2
-    mul!(L.HG, H, L.G)
+    if identityH
+        @. L.HG = L.G
+    else
+        calcHarr!(L.HG, L.G, L.cone)
+        @. L.HG *= mu
+    end
     mul!(L.GHG, L.G', L.HG)
     mul!(L.GHGQ2, L.GHG, L.Q2)
     mul!(L.Q2GHGQ2, L.Q2', L.GHGQ2)
diff --git a/src/mathoptinterface.jl b/src/mathoptinterface.jl
index d5c364293..baf3f190b 100644
--- a/src/mathoptinterface.jl
+++ b/src/mathoptinterface.jl
@@ -400,7 +400,7 @@ function MOI.optimize!(moiopt::HypatiaOptimizer)
     # TODO make it optional
     (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = preprocess_data(c, A, b, G, useQR=true)
 
-    L = QRSymmCache(c1, A1, b1, G1, h, Q2, RiQ1)
+    L = QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1)
     load_data!(opt, c1, A1, b1, G1, h, cone, L)
 
     solve!(opt)
diff --git a/src/nativeinterface.jl b/src/nativeinterface.jl
index 8231d8b1e..ad4d2a70e 100644
--- a/src/nativeinterface.jl
+++ b/src/nativeinterface.jl
@@ -316,7 +316,6 @@ function solve!(opt::Optimizer)
     loadpnt!(cone, ls_ts)
     # gradient evaluations at ls_ts of the barrier function for K
     g = similar(ts)
-    H = Matrix(1.0I, q, q)
     # helper arrays for residuals, right-hand-sides, and search directions
     tmp_tx = similar(tx)
     tmp_tx2 = similar(tx)
@@ -335,7 +334,7 @@ function solve!(opt::Optimizer)
     @. tx = 0.0
     @. ty = -b
     @. tz = -h
-    solvelinsys3!(tx, ty, tz, H, opt.L)
+    solvelinsys3!(tx, ty, tz, 1.0, opt.L, identityH=true)
     @. ts = -tz
 
     # from ts, step along interior direction of cone
@@ -501,8 +500,7 @@ function solve!(opt::Optimizer)
         # calculate prediction direction
         @. tmp_ts = tmp_tz
         @. tmp_tz = -tz
-        calcHarr!(H, Matrix(mu*I, q, q), cone) # TODO rewrite function to just return the H (or vector of H_k)
-        (tmp_kap, tmp_tau) = solvelinsys6!(tmp_tx, tmp_ty, tmp_tz, tmp_ts, -kap, kap + cx + by + hz, mu, tau, H, opt.L)
+        (tmp_kap, tmp_tau) = solvelinsys6!(tmp_tx, tmp_ty, tmp_tz, tmp_ts, -kap, kap + cx + by + hz, mu, tau, opt.L)
 
         # determine step length alpha by line search
         alpha = alphapred
@@ -590,8 +588,7 @@ function solve!(opt::Optimizer)
             calcg!(g, cone)
             @. tmp_tz = -tz - mu*g
             @. tmp_ts = 0.0
-            calcHarr!(H, Matrix(mu*I, q, q), cone) # TODO rewrite function to just return the H (or vector of H_k)
-            (tmp_kap, tmp_tau) = solvelinsys6!(tmp_tx, tmp_ty, tmp_tz, tmp_ts, -kap + mu/tau, 0.0, mu, tau, H, opt.L)
+            (tmp_kap, tmp_tau) = solvelinsys6!(tmp_tx, tmp_ty, tmp_tz, tmp_ts, -kap + mu/tau, 0.0, mu, tau, opt.L)
 
             # determine step length alpha by line search
             alpha = opt.alphacorr
diff --git a/test/native.jl b/test/native.jl
index 0a0caf201..1bc051dc2 100644
--- a/test/native.jl
+++ b/test/native.jl
@@ -177,7 +177,7 @@ function _rsoc1(verbose::Bool, lscachetype)
     cone = Hypatia.Cone([Hypatia.RotatedSecondOrderCone(4)], [1:4])
     r = fullsolve(opt, c, A, b, G, h, cone)
     @test r.status == :Optimal
-    @test r.niters <= 15
+    @test r.niters <= 20
     @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4
     @test r.pobj ≈ -sqrt(2) atol=1e-4 rtol=1e-4
     @test r.x[3:4] ≈ [1, 1]/sqrt(2) atol=1e-4 rtol=1e-4
diff --git a/test/runtests.jl b/test/runtests.jl
index 9d37254aa..352ad81e0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -14,7 +14,10 @@ using SparseArrays
 function fullsolve(opt::Hypatia.Optimizer, c, A, b, G, h, cone) # TODO handle lscachetype
     Hypatia.check_data(c, A, b, G, h, cone)
     (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true)
-    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, Q2, RiQ1)
+
+    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1)
+    # L = Hypatia.NaiveCache(c1, A1, b1, G1, h, cone)
+    
     Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L)
 
     Hypatia.solve!(opt)