diff --git a/src/algorithms/admm_one_level.jl b/src/algorithms/admm_one_level.jl
index a80fdb8..0048280 100644
--- a/src/algorithms/admm_one_level.jl
+++ b/src/algorithms/admm_one_level.jl
@@ -6,8 +6,8 @@ function admm_one_level(
     info = mod.info
     sol = mod.solution
 
-    sqrt_d = sqrt(mod.nvar)
-    OUTER_TOL = sqrt_d*(par.outer_eps) #adjusted outer loop tolerance
+    info.primtol = par.RELTOL
+    info.dualtol = par.RELTOL
 
     fill!(info, 0)
     info.mismatch = Inf
@@ -23,11 +23,11 @@ function admm_one_level(
 
     if par.verbose > 0
         admm_update_residual(env, mod, device)
-        @printf("%8s  %10s  %10s  %10s  %10s %10s %10s\n",
-                "Iter", "Objval", "Auglag", "PrimRes", "PrimTol", "DualRes", "DualTol")
+        @printf("%8s  %10s  %10s  %10s  %10s  %10s  %10s %10s %10s %10s\n",
+                "Iter", "Objval", "Auglag", "PrimScal", "PrimRes", "PrimTol", "DualScal", "DualRes", "DualTol", "rho")
 
-        @printf("%8d  %10.3e  %10.3e  %10.3e  %10.3e %10.3e  %10.3e\n",
-                info.outer, info.objval, info.auglag, info.mismatch, OUTER_TOL, info.dualres, OUTER_TOL*norm(sol.rho)/sqrt_d)
+        @printf("%8d  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e %10.3e  %10.3e  %10.3e\n",
+                info.outer, info.objval, info.auglag, info.primsca, info.primres, info.primtol, info.dualsca, info.dualres, info.dualtol, sol.rho[1])
     end
 
     info.status = :IterationLimit
@@ -51,22 +51,31 @@ function admm_one_level(
 
             if par.verbose > 0
                 if (info.cumul % 50) == 0
-                    @printf("%8s  %10s  %10s  %10s  %10s  %10s  %10s\n",
-                            "Iter", "Objval", "Auglag", "PrimRes", "PrimTol", "DualRes", "DualTol")
+                    @printf("%8s  %10s  %10s  %10s  %10s  %10s  %10s %10s %10s %10s\n",
+                            "Iter", "Objval", "Auglag", "PrimScal", "PrimRes", "PrimTol", "DualScal", "DualRes", "DualTol", "rho")
                 end
 
-                @printf("%8d  %10.3e  %10.3e  %10.3e  %10.3e %10.3e  %10.3e\n",
-                        info.outer, info.objval, info.auglag, info.mismatch, OUTER_TOL, info.dualres, OUTER_TOL*norm(sol.rho)/sqrt_d)
+                @printf("%8d  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e %10.3e  %10.3e  %10.3e\n",
+                info.outer, info.objval, info.auglag, info.primsca, info.primres, info.primtol, info.dualsca, info.dualres, info.dualtol, sol.rho[1])
             end
 
         end # while inner
 
         # mismatch: x-xbar
-        if info.mismatch <= OUTER_TOL && info.dualres <= OUTER_TOL*norm(sol.rho, device)/sqrt_d
+        if info.primres <= info.primtol # && info.dualres <= info.dualtol
             info.status = :Solved
             break
         end
 
+        # residual balancing
+        if par.rb_switch
+            if info.primres > par.rb_beta1 * info.dualres
+                sol.rho .= min.(sol.rho * par.rb_tau, 1e+16)
+            elseif par.rb_beta2 * info.primres < info.dualres
+                sol.rho .= max.(sol.rho / par.rb_tau, 1e-16)
+            end
+        end
+
     end # while outer
     end # @timed
 
diff --git a/src/algorithms/admm_two_level.jl b/src/algorithms/admm_two_level.jl
index 559ffbc..4fb9b49 100644
--- a/src/algorithms/admm_two_level.jl
+++ b/src/algorithms/admm_two_level.jl
@@ -42,7 +42,7 @@ function admm_two_level(
             admm_update_residual(env, mod, device)
 
             # an adjusting termination criteria for inner loop (i.e., inner loop is not solved to exact)
-            info.eps_pri = sqrt_d / (2500*info.outer)
+            info.eps_pri = min(sqrt_d / (2500*info.outer) / info.primsca, info.primres)
 
             if par.verbose > 0
                 if (info.cumul % 50) == 0
diff --git a/src/interface/solve_qpsub.jl b/src/interface/solve_qpsub.jl
index 97f6c4e..7759a0d 100644
--- a/src/interface/solve_qpsub.jl
+++ b/src/interface/solve_qpsub.jl
@@ -34,7 +34,8 @@ function solve_qpsub(
     use_linelimit = true,
     use_projection = false,
     tight_factor = 1.0,
-    outer_eps = 2 * 1e-4,
+    ABSTOL = 2 * 1e-4,
+    RELTOL = 1e-4,
     gpu_no = 0,
     verbose = 1,
     onelevel = true,
@@ -102,7 +103,8 @@ function solve_qpsub(
 
     env.params.scale = scale
     env.params.obj_scale = obj_scale
-    env.params.outer_eps = outer_eps
+    env.params.ABSTOL = ABSTOL
+    env.params.RELTOL = RELTOL
     env.params.outer_iterlim = outer_iterlim
     env.params.inner_iterlim = inner_iterlim
     env.params.shmem_size = sizeof(Float64) * (16 * mod.n + 4 * mod.n^2 + 178) + sizeof(Int) * (4 * mod.n)
diff --git a/src/models/acopf/acopf_admm_update_residual_cpu.jl b/src/models/acopf/acopf_admm_update_residual_cpu.jl
index 9d8a29f..0eb2e11 100644
--- a/src/models/acopf/acopf_admm_update_residual_cpu.jl
+++ b/src/models/acopf/acopf_admm_update_residual_cpu.jl
@@ -13,7 +13,8 @@ The primal and dual residuals are stored in `mod.solution.rp` and `mod.solution.
 function admm_update_residual(
     env::AdmmEnv{Float64,Array{Float64,1},Array{Int,1},Array{Float64,2}},
     mod::AbstractOPFModel{Float64,Array{Float64,1},Array{Int,1},Array{Float64,2}},
-    device::Nothing=nothing
+    device::Nothing=nothing;
+    normalized=true
 )
     sol, info, data, par, grid_data = mod.solution, mod.info, env.data, env.params, mod.grid_data
 
@@ -21,8 +22,18 @@ function admm_update_residual(
     sol.rd .= sol.z_curr .- sol.z_prev
     sol.Ax_plus_By .= sol.rp .- sol.z_curr
 
+    info.primsca = max(norm(sol.u_curr), norm(sol.v_curr), norm(sol.z_curr))
+    info.dualsca = norm(sol.l_curr)
     info.primres = norm(sol.rp)
     info.dualres = norm(sol.rd)
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.primsca
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.dualsca
+    if normalized
+        info.primres /= info.primsca
+        info.dualres /= info.dualsca
+        info.primtol /= info.primsca
+        info.dualtol /= info.dualsca
+    end
     info.norm_z_curr = norm(sol.z_curr)
     info.mismatch = norm(sol.Ax_plus_By)
     info.objval = sum(data.generators[g].coeff[data.generators[g].n-2]*(grid_data.baseMVA*sol.u_curr[mod.gen_start+2*(g-1)])^2 +
diff --git a/src/models/acopf/acopf_admm_update_residual_gpu.jl b/src/models/acopf/acopf_admm_update_residual_gpu.jl
index 6754247..d5a8808 100644
--- a/src/models/acopf/acopf_admm_update_residual_gpu.jl
+++ b/src/models/acopf/acopf_admm_update_residual_gpu.jl
@@ -11,17 +11,28 @@ end
 function admm_update_residual(
     env::AdmmEnv{Float64,CuArray{Float64,1},CuArray{Int,1},CuArray{Float64,2}},
     mod::AbstractOPFModel{Float64,CuArray{Float64,1},CuArray{Int,1},CuArray{Float64,2}},
-    device::Nothing=nothing
+    device::Nothing=nothing;
+    normalized=true
 )
-    sol, info = mod.solution, mod.info
+    sol, info, par = mod.solution, mod.info, env.params
 
     @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) compute_primal_residual_kernel(mod.nvar, sol.rp, sol.u_curr, sol.v_curr, sol.z_curr)
     @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) vector_difference(mod.nvar, sol.rd, sol.z_curr, sol.z_prev)
     CUDA.synchronize()
     CUDA.@sync @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) vector_difference(mod.nvar, sol.Ax_plus_By, sol.rp, sol.z_curr)
 
+    info.primsca = max(CUDA.norm(sol.u_curr), CUDA.norm(sol.v_curr), CUDA.norm(sol.z_curr))
+    info.dualsca = CUDA.norm(sol.l_curr)
     info.primres = CUDA.norm(sol.rp)
     info.dualres = CUDA.norm(sol.rd)
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.primsca
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.dualsca
+    if normalized
+        info.primres /= info.primsca
+        info.dualres /= info.dualsca
+        info.primtol /= info.primsca
+        info.dualtol /= info.dualsca
+    end
     info.norm_z_curr = CUDA.norm(sol.z_curr)
     info.mismatch = CUDA.norm(sol.Ax_plus_By)
 
diff --git a/src/models/acopf/acopf_admm_update_residual_ka.jl b/src/models/acopf/acopf_admm_update_residual_ka.jl
index 429c48f..51628f6 100644
--- a/src/models/acopf/acopf_admm_update_residual_ka.jl
+++ b/src/models/acopf/acopf_admm_update_residual_ka.jl
@@ -13,9 +13,10 @@ end
 function admm_update_residual(
     env::AdmmEnv,
     mod::AbstractOPFModel,
-    device
+    device;
+    normalized=true
 )
-    sol, info = mod.solution, mod.info
+    sol, info, par = mod.solution, mod.info, env.params
     nblk_nvar = div(mod.nvar-1, 64)+1
     compute_primal_residual_kernel_ka(device,64,64*nblk_nvar)(
         mod.nvar, sol.rp, sol.u_curr, sol.v_curr, sol.z_curr
@@ -30,8 +31,18 @@ function admm_update_residual(
     )
     KA.synchronize(device)
 
+    info.primsca = max(norm(sol.u_curr, device), norm(sol.v_curr, device), norm(sol.z_curr, device))
+    info.dualsca = norm(sol.l_curr, device)
     info.primres = norm(sol.rp, device)
     info.dualres = norm(sol.rd, device)
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.primsca
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.dualsca
+    if normalized
+        info.primres /= info.primsca
+        info.dualres /= info.dualsca
+        info.primtol /= info.primsca
+        info.dualtol /= info.dualsca
+    end
     info.norm_z_curr = norm(sol.z_curr, device)
     info.mismatch = norm(sol.Ax_plus_By, device)
 
diff --git a/src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl b/src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl
index d1fde6c..d109975 100644
--- a/src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl
+++ b/src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl
@@ -3,15 +3,19 @@ function admm_update_residual(
     mod::ModelMpacopf{Float64,Array{Float64,1},Array{Int,1},Array{Float64,2}},
     device::Nothing=nothing
 )
-    info = mod.info
+    info, par = mod.info, env.params
 
     info.primres = 0.0
     info.dualres = 0.0
+    info.primsca = 0.0
+    info.dualsca = 0.0
     info.norm_z_curr = 0.0
     info.mismatch = 0.0
 
     for i=1:mod.len_horizon
-        admm_update_residual(env, mod.models[i])
+        admm_update_residual(env, mod.models[i], normalized = false)
+        info.primsca = max(info.primsca, mod.models[i].info.primsca)
+        info.dualsca = max(info.dualsca, mod.models[i].info.dualsca)
     #=
         info.primres += mod.models[i].info.primres^2
         info.dualres += mod.models[i].info.dualres^2
@@ -31,6 +35,8 @@ function admm_update_residual(
         mod.models[i].info.dualres = sqrt(mod.models[i].info.dualres^2 + norm(sol_ramp.rd)^2)
         mod.models[i].info.norm_z_curr = sqrt(mod.models[i].info.norm_z_curr^2 + norm(sol_ramp.z_curr)^2)
         mod.models[i].info.mismatch = sqrt(mod.models[i].info.mismatch^2 + norm(sol_ramp.Ax_plus_By)^2)
+        info.primsca = max(info.primsca, norm(sol_ramp.u_curr), norm(v_curr), norm(sol_ramp.z_curr))
+        info.dualsca = max(info.dualsca, norm(sol_ramp.l_curr))
     #=
         info.primres += norm(sol_ramp.rp)^2
         info.dualres += norm(sol_ramp.rd)^2
@@ -45,6 +51,10 @@ function admm_update_residual(
         info.norm_z_curr = max(info.norm_z_curr, mod.models[i].info.norm_z_curr)
         info.mismatch = max(info.mismatch, mod.models[i].info.mismatch)
     end
+    info.primres /= info.primsca
+    info.dualres /= info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
 
     #=
     info.primres = sqrt(info.primres)
diff --git a/src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl b/src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl
index 51fc4c5..a0345a9 100644
--- a/src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl
+++ b/src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl
@@ -16,15 +16,19 @@ function admm_update_residual(
     mod::ModelMpacopf{Float64,CuArray{Float64,1},CuArray{Int,1},CuArray{Float64,2}},
     device::Nothing=nothing
 )
-    info = mod.info
+    info, par = mod.info, env.params
 
     info.primres = 0.0
     info.dualres = 0.0
+    info.primsca = 0.0
+    info.dualsca = 0.0
     info.norm_z_curr = 0.0
     info.mismatch = 0.0
 
     for i=1:mod.len_horizon
-        admm_update_residual(env, mod.models[i])
+        admm_update_residual(env, mod.models[i], normalized = false)
+        info.primsca = max(info.primsca, mod.models[i].info.primsca)
+        info.dualsca = max(info.dualsca, mod.models[i].info.dualsca)
     #=
         info.primres += mod.models[i].info.primres^2
         info.dualres += mod.models[i].info.dualres^2
@@ -44,6 +48,8 @@ function admm_update_residual(
         mod.models[i].info.dualres = sqrt(mod.models[i].info.dualres^2 + CUDA.norm(sol_ramp.rd)^2)
         mod.models[i].info.norm_z_curr = sqrt(mod.models[i].info.norm_z_curr^2 + CUDA.norm(sol_ramp.z_curr)^2)
         mod.models[i].info.mismatch = sqrt(mod.models[i].info.mismatch^2 + CUDA.norm(sol_ramp.Ax_plus_By)^2)
+        info.primsca = max(info.primsca, CUDA.norm(sol_ramp.u_curr), CUDA.norm(sol_ramp.v_curr), CUDA.norm(sol_ramp.z_curr))
+        info.dualsca = max(info.dualsca, CUDA.norm(sol_ramp.l_curr))
     #=
         info.primres += norm(sol_ramp.rp)^2
         info.dualres += norm(sol_ramp.rd)^2
@@ -58,6 +64,10 @@ function admm_update_residual(
         info.norm_z_curr = max(info.norm_z_curr, mod.models[i].info.norm_z_curr)
         info.mismatch = max(info.mismatch, mod.models[i].info.mismatch)
     end
+    info.primres /= info.primsca
+    info.dualres /= info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
 
     #=
     info.primres = sqrt(info.primres)
diff --git a/src/models/qpsub/qpsub_admm_update_residual_cpu.jl b/src/models/qpsub/qpsub_admm_update_residual_cpu.jl
index 607fd8b..561bafa 100644
--- a/src/models/qpsub/qpsub_admm_update_residual_cpu.jl
+++ b/src/models/qpsub/qpsub_admm_update_residual_cpu.jl
@@ -19,8 +19,12 @@ function admm_update_residual(
     sol.rd .= sol.rho .* (sol.v_curr - mod.v_prev)#from Boyd's single-level admm
     sol.Ax_plus_By .= sol.rp #x-xbar
 
-    info.primres = norm(sol.rp)
-    info.dualres = norm(sol.rd)
+    info.primsca = max(norm(sol.u_curr), norm(sol.v_curr))
+    info.dualsca = norm(sol.l_curr)
+    info.primres = norm(sol.rp) / info.primsca
+    info.dualres = norm(sol.rd) / info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
     info.mismatch = norm(sol.Ax_plus_By)
 
 
diff --git a/src/models/qpsub/qpsub_admm_update_residual_gpu.jl b/src/models/qpsub/qpsub_admm_update_residual_gpu.jl
index 98e7f33..8c27a1f 100644
--- a/src/models/qpsub/qpsub_admm_update_residual_gpu.jl
+++ b/src/models/qpsub/qpsub_admm_update_residual_gpu.jl
@@ -42,11 +42,12 @@ function admm_update_residual(
 
     @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) copy_data_kernel(mod.nvar, sol.Ax_plus_By, sol.rp) # from gpu utility
 
-
-    info.primres = CUDA.norm(sol.rp)
-
-    info.dualres = CUDA.norm(sol.rd)
-
+    info.primsca = max(CUDA.norm(sol.u_curr), CUDA.norm(sol.v_curr))
+    info.dualsca = CUDA.norm(sol.l_curr)
+    info.primres = CUDA.norm(sol.rp) / info.primsca
+    info.dualres = CUDA.norm(sol.rd) / info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
     info.mismatch = CUDA.norm(sol.Ax_plus_By)
 
     return
diff --git a/src/models/qpsub/qpsub_admm_update_residual_ka.jl b/src/models/qpsub/qpsub_admm_update_residual_ka.jl
index 082bb1a..9500dff 100644
--- a/src/models/qpsub/qpsub_admm_update_residual_ka.jl
+++ b/src/models/qpsub/qpsub_admm_update_residual_ka.jl
@@ -51,11 +51,12 @@ function admm_update_residual(
     ) # from gpu utility
     KA.synchronize(device)
 
-
-    info.primres = norm(sol.rp, device)
-
-    info.dualres = norm(sol.rd, device)
-
+    info.primsca = max(norm(sol.u_curr, device), norm(sol.v_curr, device))
+    info.dualsca = norm(sol.l_curr, device)
+    info.primres = norm(sol.rp, device) / info.primsca
+    info.dualres = norm(sol.rd, device) / info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
     info.mismatch = norm(sol.Ax_plus_By, device)
 
     return
diff --git a/src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_gpu.jl b/src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_gpu.jl
index ce5e00b..891f303 100644
--- a/src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_gpu.jl
+++ b/src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_gpu.jl
@@ -51,6 +51,8 @@ function auglag_linelimit_qpsub(Hs, l_curr, rho, u_curr, v_curr, z_curr, YffR, Y
     Atron = CuDynamicSharedArray(Float64, (6,6), (5*n + n^2 + 136)*sizeof(Float64))
     btron = CuDynamicSharedArray(Float64, 6, (5*n + n^2 + 172)*sizeof(Float64))
 
+    fill!(Ctron,0.0)
+    CUDA.sync_threads()
 
     #initialization: variable wrt branch structure wrt Exatron
     x[1] = 0.0
@@ -108,7 +110,6 @@ function auglag_linelimit_qpsub(Hs, l_curr, rho, u_curr, v_curr, z_curr, YffR, Y
     inv12 = -LH_1h[lineidx,2]/prod
     inv21 = -LH_1i[lineidx,1]/prod
     inv22 = LH_1h[lineidx,1]/prod
-    fill!(Ctron,0.0)
     Ctron[1,1] = 1.0
     Ctron[2,2] = 1.0
     Ctron[3,1] = 0.0
@@ -162,7 +163,7 @@ function auglag_linelimit_qpsub(Hs, l_curr, rho, u_curr, v_curr, z_curr, YffR, Y
 
 
     while !terminate
-         it += 1
+        it += 1
 
         eval_A_b_auglag_branch_kernel_gpu_qpsub_red(Hbr, bbr, A_aug, Atron, btron, scale,vec_1j,vec_1k,membuf,lineidx, Ctron,dtron,RH_1j,RH_1k)
 
@@ -215,7 +216,7 @@ function auglag_linelimit_qpsub(Hs, l_curr, rho, u_curr, v_curr, z_curr, YffR, Y
 
         if it >= max_auglag #maximum iteration for auglag
             if tx == 1
-                @cuprintln("max iteration reach at block = ",lineidx, "and threads = ",tx)
+                @cuprintln("max_auglag reached for line with cnorm = ", cnorm, " for line = ", lineidx, " with mu ", mu)
             end
             terminate = true
         end
diff --git a/src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_ka.jl b/src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_ka.jl
index f319290..bc50d0a 100644
--- a/src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_ka.jl
+++ b/src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_ka.jl
@@ -53,6 +53,10 @@
     Atron = @localmem Float64 (6,6)
     btron = @localmem Float64 (6,)
 
+    fill!(Ctron,0.0)
+    fill!(A_aug, 0.0)
+    fill!(Atron, 0.0)
+    CUDA.sync_threads()
 
     #initialization: variable wrt branch structure wrt Exatron
     x[1] = 0.0
@@ -112,7 +116,6 @@
     inv12 = -LH_1h[lineidx,2]/prod
     inv21 = -LH_1i[lineidx,1]/prod
     inv22 = LH_1h[lineidx,1]/prod
-    fill!(Ctron,0.0)
     Ctron[1,1] = 1.0
     Ctron[2,2] = 1.0
     Ctron[3,1] = 0.0
@@ -166,7 +169,7 @@
 
 
     while !terminate
-         it += 1
+        it += 1
 
         eval_A_b_auglag_branch_kernel_gpu_qpsub_red(
             Hbr, bbr, A_aug, Atron, btron, scale,vec_1j,vec_1k,membuf,lineidx, Ctron,dtron,RH_1j,RH_1k,
diff --git a/src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_gpu.jl b/src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_gpu.jl
index 9b7dba1..3c86d1d 100644
--- a/src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_gpu.jl
+++ b/src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_gpu.jl
@@ -72,9 +72,12 @@ end
 
 function eval_A_b_auglag_branch_kernel_gpu_qpsub_red(Hbr, bbr, A_aug, Atron, btron, scale,vec_1j,vec_1k,membuf,lineidx, Ctron,dtron,RH_1j,RH_1k)
     tx = threadIdx().x
+
+    fill!(A_aug, 0.0)
+    fill!(Atron, 0.0)
+    CUDA.sync_threads()
+
     if tx == 1
-        fill!(A_aug, 0.0)
-        fill!(Atron, 0.0)
 
         @inbounds begin
                     for i = 3:8
diff --git a/src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_ka.jl b/src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_ka.jl
index 0cecc83..586c3d9 100644
--- a/src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_ka.jl
+++ b/src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_ka.jl
@@ -72,9 +72,6 @@ end
 function eval_A_b_auglag_branch_kernel_gpu_qpsub_red(Hbr, bbr, A_aug, Atron, btron, scale,vec_1j,vec_1k,membuf,lineidx, Ctron,dtron,RH_1j,RH_1k, I, J)
     tx = J
     if tx == 1
-        fill!(A_aug, 0.0)
-        fill!(Atron, 0.0)
-
         @inbounds begin
                     for i = 3:8
                         for j = 3:8
diff --git a/src/models/qpsub/qpsub_init_solution_cpu.jl b/src/models/qpsub/qpsub_init_solution_cpu.jl
index 5b24fb0..4e0891d 100644
--- a/src/models/qpsub/qpsub_init_solution_cpu.jl
+++ b/src/models/qpsub/qpsub_init_solution_cpu.jl
@@ -8,7 +8,7 @@
 function init_solution!(
     model::ModelQpsub{Float64,Array{Float64,1},Array{Int,1},Array{Float64,2}},
     sol::Solution{Float64,Array{Float64,1}},
-    rho_pq::Float64, rho_va::Float64,
+    rho_pq::Float64, rho_va::Float64, warm_start::Bool,
     device::Nothing
     )
 
@@ -27,40 +27,44 @@ function init_solution!(
     YftR = model.grid_data.YftR; YftI = model.grid_data.YftI
     YtfR = model.grid_data.YtfR; YtfI = model.grid_data.YtfI
 
-    fill!(sol, 0.0)
-    fill!(model.lambda, 0.0)
+    if !warm_start
+        fill!(sol, 0.0)
+        fill!(model.lambda, 0.0)
 
-    #qpsub var
-    sol.rho .= rho_pq
+        #qpsub var
+        sol.rho .= rho_pq
 
-    for g=1:ngen
-        pg_idx = model.gen_start + 2*(g-1)
-        sol.v_curr[pg_idx] = 0.5*(model.qpsub_pgmin[g] + model.qpsub_pgmax[g])
-        sol.v_curr[pg_idx+1] = 0.5*(model.qpsub_qgmin[g] + model.qpsub_qgmax[g])
+        for g=1:ngen
+            pg_idx = model.gen_start + 2*(g-1)
+            sol.v_curr[pg_idx] = 0.5*(model.qpsub_pgmin[g] + model.qpsub_pgmax[g])
+            sol.v_curr[pg_idx+1] = 0.5*(model.qpsub_qgmin[g] + model.qpsub_qgmax[g])
+        end
     end
 
-    for l=1:nline
-
-        pij_idx = model.line_start + 8*(l-1)
-
-        supY = [YftR[l] YftI[l] YffR[l] 0 0 0;
-               -YftI[l] YftR[l] -YffI[l] 0 0 0;
-               YtfR[l] -YtfI[l] 0 YttR[l] 0 0;
-               -YtfI[l] -YtfR[l] 0 -YttI[l] 0 0] #wijR, wijI, wi, wj, theta_i, theta_j
 
+    for l=1:nline
         #initialize sqp_line
         sqp_line[:,l] = (ls[l,:] + us[l,:])/2  # order |w_ijR  | w_ijI |  wi(ij) | wj(ji) |  thetai(ij) |  thetaj(ji)|
 
-        sol.v_curr[pij_idx] = dot(supY[1,:],sqp_line[:,l])  #p_ij
-        sol.v_curr[pij_idx+1] = dot(supY[2,:],sqp_line[:,l]) #q_ij
-        sol.v_curr[pij_idx+2] = dot(supY[3,:],sqp_line[:,l]) #p_ji
-        sol.v_curr[pij_idx+3] = dot(supY[4,:],sqp_line[:,l]) #q_ji
-        sol.v_curr[pij_idx+4] = sqp_line[3,l] #w_i
-        sol.v_curr[pij_idx+5] = sqp_line[4,l] #w_j
-        sol.v_curr[pij_idx+6] = sqp_line[5,l] #theta_i
-        sol.v_curr[pij_idx+7] = sqp_line[6,l] #theta_j
+        if !warm_start
+            pij_idx = model.line_start + 8*(l-1)
+    
+            supY = [YftR[l] YftI[l] YffR[l] 0 0 0;
+                   -YftI[l] YftR[l] -YffI[l] 0 0 0;
+                   YtfR[l] -YtfI[l] 0 YttR[l] 0 0;
+                   -YtfI[l] -YtfR[l] 0 -YttI[l] 0 0] #wijR, wijI, wi, wj, theta_i, theta_j
+
+            sol.v_curr[pij_idx] = dot(supY[1,:],sqp_line[:,l])  #p_ij
+            sol.v_curr[pij_idx+1] = dot(supY[2,:],sqp_line[:,l]) #q_ij
+            sol.v_curr[pij_idx+2] = dot(supY[3,:],sqp_line[:,l]) #p_ji
+            sol.v_curr[pij_idx+3] = dot(supY[4,:],sqp_line[:,l]) #q_ji
+            sol.v_curr[pij_idx+4] = sqp_line[3,l] #w_i
+            sol.v_curr[pij_idx+5] = sqp_line[4,l] #w_j
+            sol.v_curr[pij_idx+6] = sqp_line[5,l] #theta_i
+            sol.v_curr[pij_idx+7] = sqp_line[6,l] #theta_j
 
-        sol.rho[pij_idx:pij_idx+7] .= rho_va
+            sol.rho[pij_idx:pij_idx+7] .= rho_va
+        end
     end
 
     return
diff --git a/src/models/qpsub/qpsub_init_solution_gpu.jl b/src/models/qpsub/qpsub_init_solution_gpu.jl
index 7f74657..5cc66aa 100644
--- a/src/models/qpsub/qpsub_init_solution_gpu.jl
+++ b/src/models/qpsub/qpsub_init_solution_gpu.jl
@@ -26,7 +26,7 @@ function init_branch_bus_kernel_qpsub(n::Int, line_start::Int, rho_va::Float64,
     YtfR::CuDeviceArray{Float64,1}, YtfI::CuDeviceArray{Float64,1},
     YttR::CuDeviceArray{Float64,1}, YttI::CuDeviceArray{Float64,1},
     us::CuDeviceArray{Float64,2}, ls::CuDeviceArray{Float64,2}, sqp_line::CuDeviceArray{Float64,2},
-    v::CuDeviceArray{Float64,1}, rho::CuDeviceArray{Float64,1}, supY::CuDeviceArray{Float64,2}
+    v::CuDeviceArray{Float64,1}, rho::CuDeviceArray{Float64,1}, supY::CuDeviceArray{Float64,2}, warm_start::Bool
 )
     l = threadIdx().x + (blockDim().x * (blockIdx().x - 1))
     if l <= n
@@ -53,16 +53,18 @@ function init_branch_bus_kernel_qpsub(n::Int, line_start::Int, rho_va::Float64,
         supY[4*(l-1) + 4,6] = -YttI[l]
 
 
-        v[pij_idx] = YftR[l]*sqp_line[1,l] + YftI[l]*sqp_line[2,l] + YffR[l]*sqp_line[3,l]  #CUBLAS.dot(4, supY[1,:],sqp_line[:,l])  #p_ij
-        v[pij_idx+1] = -YftI[l]*sqp_line[1,l] + YftR[l]*sqp_line[2,l] - YffI[l]*sqp_line[3,l] #CUBLAS.dot(4, supY[2,:],sqp_line[:,l]) #q_ij
-        v[pij_idx+2] = YtfR[l]*sqp_line[1,l] - YtfI[l]*sqp_line[2,l] + YttR[l]*sqp_line[4,l] #CUBLAS.dot(4, supY[3,:],sqp_line[:,l]) #p_ji
-        v[pij_idx+3] = -YtfI[l]*sqp_line[1,l] - YtfR[l]*sqp_line[2,l] - YttI[l]*sqp_line[4,l] #CUBLAS.dot(4, supY[4,:],sqp_line[:,l]) #q_ji
-        v[pij_idx+4] = sqp_line[3,l] #w_i
-        v[pij_idx+5] = sqp_line[4,l] #w_j
-        v[pij_idx+6] = sqp_line[5,l] #theta_i
-        v[pij_idx+7] = sqp_line[6,l] #theta_j
+        if !warm_start
+            v[pij_idx] = YftR[l]*sqp_line[1,l] + YftI[l]*sqp_line[2,l] + YffR[l]*sqp_line[3,l]  #CUBLAS.dot(4, supY[1,:],sqp_line[:,l])  #p_ij
+            v[pij_idx+1] = -YftI[l]*sqp_line[1,l] + YftR[l]*sqp_line[2,l] - YffI[l]*sqp_line[3,l] #CUBLAS.dot(4, supY[2,:],sqp_line[:,l]) #q_ij
+            v[pij_idx+2] = YtfR[l]*sqp_line[1,l] - YtfI[l]*sqp_line[2,l] + YttR[l]*sqp_line[4,l] #CUBLAS.dot(4, supY[3,:],sqp_line[:,l]) #p_ji
+            v[pij_idx+3] = -YtfI[l]*sqp_line[1,l] - YtfR[l]*sqp_line[2,l] - YttI[l]*sqp_line[4,l] #CUBLAS.dot(4, supY[4,:],sqp_line[:,l]) #q_ji
+            v[pij_idx+4] = sqp_line[3,l] #w_i
+            v[pij_idx+5] = sqp_line[4,l] #w_j
+            v[pij_idx+6] = sqp_line[5,l] #theta_i
+            v[pij_idx+7] = sqp_line[6,l] #theta_j
 
-        rho[pij_idx:pij_idx+7] .= rho_va
+            rho[pij_idx:pij_idx+7] .= rho_va
+        end
 
 
     end
@@ -73,25 +75,24 @@ end
 function init_solution!(
     model::ModelQpsub{Float64,CuArray{Float64,1},CuArray{Int,1},CuArray{Float64,2}},
     sol::Solution{Float64,CuArray{Float64,1}},
-    rho_pq::Float64, rho_va::Float64,
-    device::Nothing
+    rho_pq::Float64, rho_va::Float64, warm_start::Bool,
+    device::Nothing; 
     )
 
     data = model.grid_data
 
-    fill!(sol, 0.0)
-    fill!(model.lambda, 0.0)
-
-    #qpsub var
-    sol.rho .= rho_pq
-
-
-    @cuda threads=64 blocks=(div(data.ngen-1,64)+1) init_generator_kernel_qpsub(data.ngen, model.gen_start,
-                    model.qpsub_pgmax, model.qpsub_pgmin, model.qpsub_qgmax, model.qpsub_qgmin, sol.v_curr)
-
+    if !warm_start
+        fill!(sol, 0.0)
+        fill!(model.lambda, 0.0)
+    
+        #qpsub var
+        sol.rho .= rho_pq
+        @cuda threads=64 blocks=(div(data.ngen-1,64)+1) init_generator_kernel_qpsub(data.ngen, model.gen_start,
+                        model.qpsub_pgmax, model.qpsub_pgmin, model.qpsub_qgmax, model.qpsub_qgmin, sol.v_curr)
+    end
     @cuda threads=64 blocks=(div(data.nline-1,64)+1) init_branch_bus_kernel_qpsub(data.nline, model.line_start, rho_va,
                     data.YffR, data.YffI, data.YftR, data.YftI,
-                    data.YtfR, data.YtfI, data.YttR, data.YttI, model.us, model.ls, model.sqp_line, sol.v_curr, sol.rho, model.supY)
+                    data.YtfR, data.YtfI, data.YttR, data.YttI, model.us, model.ls, model.sqp_line, sol.v_curr, sol.rho, model.supY, warm_start)
     CUDA.synchronize()
 
     return
diff --git a/src/models/qpsub/qpsub_init_solution_ka.jl b/src/models/qpsub/qpsub_init_solution_ka.jl
index b94353f..189d32b 100644
--- a/src/models/qpsub/qpsub_init_solution_ka.jl
+++ b/src/models/qpsub/qpsub_init_solution_ka.jl
@@ -26,7 +26,7 @@ end
     YtfR, YtfI,
     YttR, YttI,
     us, ls, sqp_line,
-    v, rho, supY
+    v, rho, supY, warm_start::Bool
 )
     I = @index(Group, Linear)
     J = @index(Local, Linear)
@@ -55,17 +55,18 @@ end
         supY[4*(l-1) + 4,6] = -YttI[l]
 
 
-        v[pij_idx] = YftR[l]*sqp_line[1,l] + YftI[l]*sqp_line[2,l] + YffR[l]*sqp_line[3,l]  #CUBLAS.dot(4, supY[1,:],sqp_line[:,l])  #p_ij
-        v[pij_idx+1] = -YftI[l]*sqp_line[1,l] + YftR[l]*sqp_line[2,l] - YffI[l]*sqp_line[3,l] #CUBLAS.dot(4, supY[2,:],sqp_line[:,l]) #q_ij
-        v[pij_idx+2] = YtfR[l]*sqp_line[1,l] - YtfI[l]*sqp_line[2,l] + YttR[l]*sqp_line[4,l] #CUBLAS.dot(4, supY[3,:],sqp_line[:,l]) #p_ji
-        v[pij_idx+3] = -YtfI[l]*sqp_line[1,l] - YtfR[l]*sqp_line[2,l] - YttI[l]*sqp_line[4,l] #CUBLAS.dot(4, supY[4,:],sqp_line[:,l]) #q_ji
-        v[pij_idx+4] = sqp_line[3,l] #w_i
-        v[pij_idx+5] = sqp_line[4,l] #w_j
-        v[pij_idx+6] = sqp_line[5,l] #theta_i
-        v[pij_idx+7] = sqp_line[6,l] #theta_j
-
-        rho[pij_idx:pij_idx+7] .= rho_va
+        if !warm_start
+            v[pij_idx] = YftR[l]*sqp_line[1,l] + YftI[l]*sqp_line[2,l] + YffR[l]*sqp_line[3,l]  #CUBLAS.dot(4, supY[1,:],sqp_line[:,l])  #p_ij
+            v[pij_idx+1] = -YftI[l]*sqp_line[1,l] + YftR[l]*sqp_line[2,l] - YffI[l]*sqp_line[3,l] #CUBLAS.dot(4, supY[2,:],sqp_line[:,l]) #q_ij
+            v[pij_idx+2] = YtfR[l]*sqp_line[1,l] - YtfI[l]*sqp_line[2,l] + YttR[l]*sqp_line[4,l] #CUBLAS.dot(4, supY[3,:],sqp_line[:,l]) #p_ji
+            v[pij_idx+3] = -YtfI[l]*sqp_line[1,l] - YtfR[l]*sqp_line[2,l] - YttI[l]*sqp_line[4,l] #CUBLAS.dot(4, supY[4,:],sqp_line[:,l]) #q_ji
+            v[pij_idx+4] = sqp_line[3,l] #w_i
+            v[pij_idx+5] = sqp_line[4,l] #w_j
+            v[pij_idx+6] = sqp_line[5,l] #theta_i
+            v[pij_idx+7] = sqp_line[6,l] #theta_j
 
+            rho[pij_idx:pij_idx+7] .= rho_va
+        end
 
     end
 end
@@ -73,27 +74,29 @@ end
 function init_solution!(
     model::ModelQpsub,
     sol::Solution,
-    rho_pq::Float64, rho_va::Float64,
+    rho_pq::Float64, rho_va::Float64, warm_start::Bool,
     device
     )
 
     data = model.grid_data
 
-    fill!(sol, 0.0)
-    fill!(model.lambda, 0.0)
+    if !warm_start
+        fill!(sol, 0.0)
+        fill!(model.lambda, 0.0)
 
-    #qpsub var
-    sol.rho .= rho_pq
+        #qpsub var
+        sol.rho .= rho_pq
 
 
-    nblk_gen = div(data.ngen-1,64)+1
-    ev = init_generator_kernel_qpsub_ka(device, 64, 64*nblk_gen)(data.ngen, model.gen_start,
-                    model.qpsub_pgmax, model.qpsub_pgmin, model.qpsub_qgmax, model.qpsub_qgmin, sol.v_curr)
+        nblk_gen = div(data.ngen-1,64)+1
+        ev = init_generator_kernel_qpsub_ka(device, 64, 64*nblk_gen)(data.ngen, model.gen_start,
+                        model.qpsub_pgmax, model.qpsub_pgmin, model.qpsub_qgmax, model.qpsub_qgmin, sol.v_curr)
+    end
 
     nblk_line = div(data.nline-1,64)+1
     ev = init_branch_bus_kernel_qpsub_ka(device, 64, 64*nblk_line)(data.nline, model.line_start, rho_va,
                     data.YffR, data.YffI, data.YftR, data.YftI,
-                    data.YtfR, data.YtfI, data.YttR, data.YttI, model.us, model.ls, model.sqp_line, sol.v_curr, sol.rho, model.supY)
+                    data.YtfR, data.YtfI, data.YttR, data.YttI, model.us, model.ls, model.sqp_line, sol.v_curr, sol.rho, model.supY, warm_start)
     KA.synchronize(device)
 
     return
diff --git a/src/utils/environment.jl b/src/utils/environment.jl
index 040fde6..638ec44 100644
--- a/src/utils/environment.jl
+++ b/src/utils/environment.jl
@@ -17,6 +17,10 @@ mutable struct Parameters
     rt_inc::Float64     #? not used
     rt_dec::Float64     #? not used
     eta::Float64        #? not used
+    rb_switch::Bool     # residual balancing on/off
+    rb_tau::Float64     # residual balancing multiplier
+    rb_beta1::Float64   # residual balancing criterion
+    rb_beta2::Float64   # residual balancing criterion
     verbose::Int
 
     # MPI implementation
@@ -51,6 +55,10 @@ mutable struct Parameters
         par.rt_inc = 2.0
         par.rt_dec = 2.0
         par.eta = 0.99
+        par.rb_switch = false
+        par.rb_tau = 2
+        par.rb_beta1 = 10
+        par.rb_beta2 = 10
         par.max_auglag = 50
         par.ABSTOL = 1e-6
         par.RELTOL = 1e-5
@@ -333,6 +341,10 @@ mutable struct IterationInformation{U}
     objval::Float64
     primres::Float64
     dualres::Float64
+    primtol::Float64 # normalized primal tolerance; this may be redundant to eps_pri below.
+    dualtol::Float64 # normalized dual tolerance
+    primsca::Float64 # primal normalization scalar
+    dualsca::Float64 # primal normalization scalar
     mismatch::Float64
     auglag::Float64
     eps_pri::Float64
@@ -365,6 +377,10 @@ function Base.fill!(info::IterationInformation, val)
     info.objval = val
     info.primres = val
     info.dualres = val
+    info.primtol = val
+    info.dualtol = val
+    info.primsca = val
+    info.dualsca = val
     info.mismatch = val
     info.auglag = val
     info.eps_pri = val
@@ -389,6 +405,10 @@ function Base.copy(ref::IterationInformation{ComponentInformation})
     info.objval = ref.objval
     info.primres = ref.primres
     info.dualres = ref.dualres
+    info.primtol = ref.primtol
+    info.dualtol = ref.dualtol
+    info.primsca = ref.primsca
+    info.dualsca = ref.dualsca
     info.mismatch = ref.mismatch
     info.auglag = ref.auglag
     info.eps_pri = ref.eps_pri
diff --git a/test/algorithms/qpsub_update_cpu.jl b/test/algorithms/qpsub_update_cpu.jl
index f01acca..026a694 100644
--- a/test/algorithms/qpsub_update_cpu.jl
+++ b/test/algorithms/qpsub_update_cpu.jl
@@ -229,11 +229,11 @@ end #@testset
     env3, mod3 = ExaAdmm.solve_qpsub(case, mod1.Hs, mod1.LH_1h, mod1.RH_1h,
     mod1.LH_1i, mod1.RH_1i, mod1.LH_1j, mod1.RH_1j, mod1.LH_1k, mod1.RH_1k, mod1.ls, mod1.us, mod1.qpsub_pgmax, mod1.qpsub_pgmin, mod1.qpsub_qgmax, mod1.qpsub_qgmin, mod1.qpsub_c1, mod1.qpsub_c2, mod1.qpsub_Pd, mod1.qpsub_Qd,
     initial_beta; 
-        outer_iterlim=10000, inner_iterlim=1, scale = 1e-4, obj_scale = 1, rho_pq = 4000.0, rho_va = 4000.0, verbose=0, outer_eps=2*1e-6, onelevel = true)
+        outer_iterlim=10000, inner_iterlim=1, scale = 1e-4, obj_scale = 1, rho_pq = 4000.0, rho_va = 4000.0, verbose=verbose, ABSTOL=2e-6, RELTOL=5e-5, onelevel = true)
 
   
     @test mod3.info.status == :Solved
-    @test mod3.info.outer == 5107
-    @test mod3.info.cumul == 5107
-    @test isapprox(mod3.info.objval, -21.92744641968529; atol=1e-6)     
+    @test mod3.info.outer == 5481
+    @test mod3.info.cumul == 5481
+    @test isapprox(mod3.info.objval, -21.9278463925427; atol=2e-6)
 end 
\ No newline at end of file
diff --git a/test/algorithms/qpsub_update_gpu.jl b/test/algorithms/qpsub_update_gpu.jl
index 49a7ae4..205304a 100644
--- a/test/algorithms/qpsub_update_gpu.jl
+++ b/test/algorithms/qpsub_update_gpu.jl
@@ -1,19 +1,21 @@
-@testset "Testing [x,xbar,l,res] updates and a solve for ACOPF" begin
-
 rho_pq = 20.0 #for two level 
 rho_va = 20.0 #for two level
 initial_beta = 100000.0 #for two level
 atol = 2e-6
+rtol = 5e-5
 verbose = 0
 use_gpu = true
 
 case = joinpath(INSTANCES_DIR, "case9.m")
 
 # Initialize an cpu qpsub model with default options as shell for qpsub.
-    T = Float64; TD = Array{Float64,1}; TI = Array{Int,1}; TM = Array{Float64,2}
+T = Float64; TD = Array{Float64,1}; TI = Array{Int,1}; TM = Array{Float64,2}
+
+env1 = ExaAdmm.AdmmEnv{T,TD,TI,TM}(case, rho_pq, rho_va; verbose=verbose)
+mod1 = ExaAdmm.ModelQpsub{T,TD,TI,TM}(env1)
+
+@testset "Testing [x,xbar,l,res] updates and a solve for QPsub" begin
 
-    env1 = ExaAdmm.AdmmEnv{T,TD,TI,TM}(case, rho_pq, rho_va; verbose=verbose)
-    mod1 = ExaAdmm.ModelQpsub{T,TD,TI,TM}(env1)
     sol = mod1.solution
     par = env1.params
     data = mod1.grid_data
@@ -227,132 +229,31 @@ end
     env3, mod3 = ExaAdmm.solve_qpsub(case, mod1.Hs, mod1.LH_1h, mod1.RH_1h,
     mod1.LH_1i, mod1.RH_1i, mod1.LH_1j, mod1.RH_1j, mod1.LH_1k, mod1.RH_1k, mod1.ls, mod1.us, mod1.qpsub_pgmax, mod1.qpsub_pgmin, mod1.qpsub_qgmax, mod1.qpsub_qgmin, mod1.qpsub_c1, mod1.qpsub_c2, mod1.qpsub_Pd, mod1.qpsub_Qd,
     initial_beta; 
-        outer_iterlim=10000, inner_iterlim=1, scale = 1e-4, obj_scale = 1, rho_pq = 4000.0, rho_va = 4000.0, verbose=0, outer_eps=2*1e-6, onelevel = true, use_gpu = true)
+        outer_iterlim=10000, inner_iterlim=1, scale = 1e-4, obj_scale = 1, rho_pq = 4000.0, rho_va = 4000.0, verbose=verbose, ABSTOL=atol, RELTOL=rtol, onelevel = true, use_gpu = true)
 
   
     @test mod3.info.status == :Solved
-    @test mod3.info.outer == 5107
-    @test mod3.info.cumul == 5107
-    @test isapprox(mod3.info.objval, -21.92744641968529; atol=1e-6)    
+    @test mod3.info.outer == 5481
+    @test mod3.info.cumul == 5481
+    @test isapprox(mod3.info.objval, -21.927846392542; atol=atol)    
     
     sol3 = mod3.solution
-    atol = 1e-6
-
-    dpg_sol_cpu = [-0.11550455128840183
-    0.06549995591228526
-    0.053676513010466456]
-
-    dqg_sol_cpu = [-0.007903958357205233
-    0.019785362536583685
-    0.011746408694041054]
-
-    dline_var_cpu = [-0.00630369   -0.00556137  -0.00333908   0.000404039   9.04974e-6   -8.98924e-6    0.000661247  -0.00373931  -0.00567469
-    -0.00665306   -0.00509225  -0.00938981   0.00314544   -0.000313391  -0.000223849  -0.00409375    0.00970327   0.00512255
-    -0.00675896   -0.0063758   -0.00500354   0.00109238   -1.06479e-8   -1.48199e-6   -1.06804e-8   -1.06804e-8  -0.00534789
-    -0.00637581   -0.00500357  -1.06479e-8  -1.06479e-8   -1.48608e-6   -1.06804e-8    0.00189783   -0.00534788  -0.0063758
-    -3.25492e-19   0.00531377   0.00949253   0.0201868     0.0175742     0.0178357     0.0180226     0.0180226    0.00951113
-     0.00531381    0.00949261   0.0175743    0.0175742     0.0178358     0.0180226     0.021366      0.00951121   0.00531385]
-
-     dline_fl_cpu = [ -0.115505    -0.0551048   -0.0546064    0.0536765   -0.00307677   -0.00305414   -0.0655      0.0624192    0.0598887
-     -0.00790396   0.00183362   0.00363203   0.0117464    0.000273346   0.000464936  -0.0105801   0.0108191   -0.00273044
-      0.115505     0.0546064    0.0567533   -0.0536765    0.00305414    0.0030808     0.0655     -0.0598887   -0.0603997
-     -0.00125214  -0.00363203   0.00662172  -0.00689505  -0.000464924  -0.000239002   0.0197854   0.00273044  -0.000581486]
-
-     dtheta_sol_cpu = [-3.25491739361945e-19
-     0.021365998519708898
-     0.020186769395364127
-     0.005313812154732135
-     0.009492570006911722
-     0.017574220141571076
-     0.017835751182986154
-     0.01802260774972062
-     0.009511171536802466]
-
-     dw_sol_cpu = [-0.006758956828784941
-     0.0018978329222831053
-     0.0010923793242977361
-    -0.006375802955645012
-    -0.005003557142581212
-    -1.0647877468628053e-8
-    -1.484035364757708e-6
-    -1.068039123808262e-8
-    -0.0053478835040742636]
-
-    dual_infeas_1_cpu = [-254.11001283448402
-    111.34992505088496
-    131.50745687564282
-      0.31062229996245705
-     -8.03424172781389
-     -0.3078548536034032
-     -0.07573915941618026
-      0.012603427412195068
-      0.008651826283358158
-      0.03279612499623894] 
-
-    dual_infeas_2_cpu = [-2.337593793742112
-    -0.2062996739647522
-     0.10994572761041102
-     0.031770554831216935
-     0.02745497289931917
-    -0.11873053563714611
-    -1.2564033339874354
-    -0.12073294817351554
-     0.2121498763329952
-     0.05086202591219009]
-
-    dual_infeas_3_cpu = [0.057405134084524315
-    -0.16737701814667477
-     3.670092553395319
-     0.4052710265115234
-    -0.2342336895471848
-     0.07724866682808983
-     0.07379529136930803
-     0.004080775509334406
-    -0.12229701235980252
-    -0.0018175498192630844]
-
-    dual_infeas_4_cpu = [-0.0019510899995607155
-    0.07059765612753019
-    0.07104222896868044
-   -0.0061070561232362
-   -0.17080362240460334
-    0.0027843138173676155
-    0.003450443341928419
-    0.07157328724238488
-    0.07186014778306046
-   -0.28991023775276664]
-
-    dual_infeas_5_cpu = [-4.200422692029046
-    -0.33666635407314044
-     0.6407231013041463
-     0.0753794333813387
-     0.08217499064640721
-    -0.18387848500073856
-     1.4609510208883514
-     0.26889163927860954
-    -0.13184568446925882
-     0.065799538338272
-     0.04433557788772256
-     0.06271507036787678
-     2.807937471979002
-     0.0728632119983576
-    -0.20915786828282304
-     0.03875039932653859
-     0.020549534704044203]
-
-     dual_infeas_cpu = vcat(dual_infeas_1_cpu,dual_infeas_2_cpu,dual_infeas_3_cpu,dual_infeas_4_cpu,dual_infeas_5_cpu)
-
-     lambda_cpu = [0.488894     435.795        282.155     -0.0489912    248.787        348.203         0.135763     267.164        315.305
-     0.790211       0.88808        0.816501  -0.0653646      0.576725       0.456339      0.0687841      0.210784       0.140033
-    -8.67362e-19   -8.67362e-19   -0.0       -0.0           -4.33681e-19   -1.73472e-18  -0.0           -8.67362e-19   -8.67362e-19
-    -2.60209e-18   -0.0           -0.0       -1.73472e-18   -8.67362e-19   -8.67362e-19  -8.67362e-19   -0.0           -8.67362e-19]
 
-     @test norm(dpg_sol_cpu - Array(mod3.dpg_sol),  Inf) <= atol 
-     @test norm(dqg_sol_cpu - Array(mod3.dqg_sol),  Inf) <= atol 
-     @test norm(dline_var_cpu - Array(mod3.dline_var),  Inf) <= atol 
-     @test norm(dline_fl_cpu - Array(mod3.dline_fl),  Inf) <= atol 
-     @test norm(dtheta_sol_cpu - Array(mod3.dtheta_sol),  Inf) <= atol 
-     @test norm(dw_sol_cpu - Array(mod3.dw_sol),  Inf) <= atol 
-     @test norm(dual_infeas_cpu - Array(mod3.dual_infeas),  Inf) <= atol
-     @test norm((lambda_cpu - Array(mod3.lambda))./max.(abs.(lambda_cpu),1), Inf) <= 2e-6 
+    dpg_sol_cpu = [-0.11550471288823375, 0.06549963131056268, 0.05367626157137195]
+    dqg_sol_cpu = [-0.005718588343597756, 0.01851317137973553, 0.01058363347077568]
+    dline_var_cpu = [-0.005569819650835843 -0.005054193193183796 -0.0031356731585941602 0.0003340424001857633 9.018007914709618e-6 -8.955231272232481e-6 0.0005809649468087098 -0.0035327541819175487 -0.0051673927876596205; -0.006653071419384663 -0.005109374028865239 -0.009433791106253716 0.003145428964346758 -0.0003118874694407943 -0.00022277448888749965 -0.004093726855025406 0.009746743401438716 0.005135828873903401; -0.0058992105835546605 -0.005767638775672478 -0.004596762982732082 0.0009542434497012447 -1.0647877468628053e-8 -1.492916208770773e-6 -1.068039123808262e-8 -1.068039123808262e-8 -0.004939127111986979; -0.005767628302723796 -0.004596740386993611 -1.0647877468628053e-8 -1.0647877468628053e-8 -1.4557418980203987e-6 -1.068039123808262e-8 0.0017380383383034281 -0.004939090417984399 -0.005767609350738688; -3.25491739361945e-19 0.005338385069326234 0.009542093395593967 0.0202615747235613 0.01764627076409719 0.01790665022472073 0.018092699885635936 0.018092713807391705 0.009562285691049546; 0.005338374345348201 0.009542068575890466 0.017646245194261424 0.017646257498616476 0.01790662375090709 0.018092686631961837 0.021440809620779835 0.009562257967646896 0.005338358797672126]
+    dline_fl_cpu = [-0.11550471214209485 -0.0550886536682177 -0.0545914661320966 0.05367626219021771 -0.003062020230877622 -0.0030394927033266204 -0.0654996296804065 0.06243360473855911 0.059908274544507235; -0.0057185925819239125 0.002880225309702243 0.004752099132814017 0.010583635657260775 0.00027191961424352743 0.0004625834867261366 -0.009295610035199166 0.009533344180015505 -0.003927910576173748; 0.11550471214209485 0.05459146673835657 0.056738283619246 -0.05367626219021771 0.003039493246316152 0.0030660246149609373 0.0654996296804065 -0.05990827415794294 -0.06041605732594946; -0.0034341779841658483 -0.004752098035100081 0.00542864572694711 -0.005700563960123412 -0.0004625831869150321 -0.00023773057022993016 0.018513174263915494 0.00392790777176107 0.0005539497426347278]
+    dtheta_sol_cpu = [-3.25491739361945e-19, 0.021440809620779835, 0.0202615747235613, 0.005338372737448853, 0.009542080985742217, 0.017646257818991696, 0.01790663698781391, 0.018092700108329824, 0.009562271829348221]
+    dw_sol_cpu = [-0.0058992105835546605, 0.0017380383383034281, 0.0009542434497012447, -0.005767625476378321, -0.004596751684862847, -1.0647877468628053e-8, -1.4743290533955858e-6, -1.068039123808262e-8, -0.004939108764985689]
+    dual_infeas_cpu = [-254.11036835411423, 111.34937322795656, 131.50684084986128, 0.3119989408890842, -8.034251796583263, -0.22792671306485607, -0.14821201491131938, 0.011898057499073345, 0.009455439882319458, 0.03537463279074429, -2.345388987728991, -0.18156150826834788, 0.08875478602088671, 0.03139112597637954, 0.02813068860448726, -0.11748530455855485, -1.2620894362848545, -0.10695004291180904, 0.19961093972437144, 0.05088261867417226, 0.057870735685249315, -0.1686020205956723, 3.6700753782985167, 0.3650331971871314, -0.1936042960743483, 0.0776146861999613, 0.0740166426887498, 0.004061310842689346, -0.1217101785537621, -0.0018110619226406265, -0.0019394705261958892, 0.0708845572940984, 0.07132702076591872, -0.006077871401322642, -0.169983866840654, 0.002767546200533025, 0.0034373043082680165, 0.0718559407428788, 0.07214140668385148, -0.29046051526063166, -4.20040217042528, -0.2957215255096348, 0.5993737191181183, 0.07574586733811475, 0.08238817068754835, -0.1824165829874273, 1.4673408402946906, 0.25436253640211093, -0.11631341488058167, 0.0658499676530561, 0.04476991944709829, 0.06513389823348271, 2.8151683733740307, 0.04757957575153994, -0.18020940646433298, 0.038415832071539036, 0.02118674588334765]
+    lambda_cpu = [0.4390310148492841 435.7732781204522 282.1544242381189 -0.020994217416346203 248.8021373471199 348.2250378203491 0.16468763831157662 267.1641217085404 315.2820765395896; 0.7925836307556712 0.8917093439553776 0.8196728616326903 -0.06555427114780897 0.5792170456138986 0.45857551396302043 0.06884954532846001 0.21271856587119875 0.14205944713641747; -1.734723475976807e-18 -8.673617379884035e-19 -2.168404344971009e-19 -0.0 -4.336808689942018e-19 -0.0 -8.673617379884035e-19 -8.673617379884035e-19 -8.673617379884035e-19; -1.734723475976807e-18 -0.0 -0.0 -0.0 -4.336808689942018e-19 -0.0 -1.734723475976807e-18 -0.0 -0.0]
+
+    @test norm(dpg_sol_cpu - Array(mod3.dpg_sol), Inf) <= atol
+    @test norm(dqg_sol_cpu - Array(mod3.dqg_sol), Inf) <= atol
+    @test norm(dline_var_cpu - Array(mod3.dline_var), Inf) <= atol
+    @test norm(dline_fl_cpu - Array(mod3.dline_fl), Inf) <= atol
+    @test norm(dtheta_sol_cpu - Array(mod3.dtheta_sol), Inf) <= atol
+    @test norm(dw_sol_cpu - Array(mod3.dw_sol), Inf) <= atol
+    @test norm(dual_infeas_cpu - Array(mod3.dual_infeas), Inf) <= atol
+    @test norm(lambda_cpu - Array(mod3.lambda), Inf) <= atol
 end
\ No newline at end of file
diff --git a/test/algorithms/qpsub_update_ka.jl b/test/algorithms/qpsub_update_ka.jl
index 022770a..4926119 100644
--- a/test/algorithms/qpsub_update_ka.jl
+++ b/test/algorithms/qpsub_update_ka.jl
@@ -19,6 +19,7 @@ rho_pq = 20.0 #for two level
 rho_va = 20.0 #for two level
 initial_beta = 100000.0 #for two level
 atol = 2e-6
+rtol = 5e-5
 verbose = 0
 use_gpu = true
 
@@ -252,132 +253,31 @@ end
     env3, mod3 = ExaAdmm.solve_qpsub(case, mod1.Hs, mod1.LH_1h, mod1.RH_1h,
     mod1.LH_1i, mod1.RH_1i, mod1.LH_1j, mod1.RH_1j, mod1.LH_1k, mod1.RH_1k, mod1.ls, mod1.us, mod1.qpsub_pgmax, mod1.qpsub_pgmin, mod1.qpsub_qgmax, mod1.qpsub_qgmin, mod1.qpsub_c1, mod1.qpsub_c2, mod1.qpsub_Pd, mod1.qpsub_Qd,
     initial_beta;
-        outer_iterlim=10000, inner_iterlim=1, scale = 1e-4, obj_scale = 1, rho_pq = 4000.0, rho_va = 4000.0, verbose=0, outer_eps=2*1e-6, onelevel = true, use_gpu = true, ka_device=_device)
-
+        outer_iterlim=10000, inner_iterlim=1, scale = 1e-4, obj_scale = 1, rho_pq = 4000.0, rho_va = 4000.0, verbose=verbose, ABSTOL=atol, RELTOL=rtol, onelevel = true, use_gpu = true, ka_device=_device)
 
+  
     @test mod3.info.status == :Solved
-    @test mod3.info.outer == 5107
-    @test mod3.info.cumul == 5107
-    @test isapprox(mod3.info.objval, -21.92744641968529; atol=1e-6)
-
+    @test mod3.info.outer == 5481
+    @test mod3.info.cumul == 5481
+    @test isapprox(mod3.info.objval, -21.927846392542; atol=atol)    
+    
     sol3 = mod3.solution
-    atol = 1e-6
-
-    dpg_sol_cpu = [-0.11550455128840183
-    0.06549995591228526
-    0.053676513010466456]
-
-    dqg_sol_cpu = [-0.007903958357205233
-    0.019785362536583685
-    0.011746408694041054]
-
-    dline_var_cpu = [-0.00630369   -0.00556137  -0.00333908   0.000404039   9.04974e-6   -8.98924e-6    0.000661247  -0.00373931  -0.00567469
-    -0.00665306   -0.00509225  -0.00938981   0.00314544   -0.000313391  -0.000223849  -0.00409375    0.00970327   0.00512255
-    -0.00675896   -0.0063758   -0.00500354   0.00109238   -1.06479e-8   -1.48199e-6   -1.06804e-8   -1.06804e-8  -0.00534789
-    -0.00637581   -0.00500357  -1.06479e-8  -1.06479e-8   -1.48608e-6   -1.06804e-8    0.00189783   -0.00534788  -0.0063758
-    -3.25492e-19   0.00531377   0.00949253   0.0201868     0.0175742     0.0178357     0.0180226     0.0180226    0.00951113
-     0.00531381    0.00949261   0.0175743    0.0175742     0.0178358     0.0180226     0.021366      0.00951121   0.00531385]
-
-     dline_fl_cpu = [ -0.115505    -0.0551048   -0.0546064    0.0536765   -0.00307677   -0.00305414   -0.0655      0.0624192    0.0598887
-     -0.00790396   0.00183362   0.00363203   0.0117464    0.000273346   0.000464936  -0.0105801   0.0108191   -0.00273044
-      0.115505     0.0546064    0.0567533   -0.0536765    0.00305414    0.0030808     0.0655     -0.0598887   -0.0603997
-     -0.00125214  -0.00363203   0.00662172  -0.00689505  -0.000464924  -0.000239002   0.0197854   0.00273044  -0.000581486]
-
-     dtheta_sol_cpu = [-3.25491739361945e-19
-     0.021365998519708898
-     0.020186769395364127
-     0.005313812154732135
-     0.009492570006911722
-     0.017574220141571076
-     0.017835751182986154
-     0.01802260774972062
-     0.009511171536802466]
-
-     dw_sol_cpu = [-0.006758956828784941
-     0.0018978329222831053
-     0.0010923793242977361
-    -0.006375802955645012
-    -0.005003557142581212
-    -1.0647877468628053e-8
-    -1.484035364757708e-6
-    -1.068039123808262e-8
-    -0.0053478835040742636]
-
-    dual_infeas_1_cpu = [-254.11001283448402
-    111.34992505088496
-    131.50745687564282
-      0.31062229996245705
-     -8.03424172781389
-     -0.3078548536034032
-     -0.07573915941618026
-      0.012603427412195068
-      0.008651826283358158
-      0.03279612499623894]
-
-    dual_infeas_2_cpu = [-2.337593793742112
-    -0.2062996739647522
-     0.10994572761041102
-     0.031770554831216935
-     0.02745497289931917
-    -0.11873053563714611
-    -1.2564033339874354
-    -0.12073294817351554
-     0.2121498763329952
-     0.05086202591219009]
-
-    dual_infeas_3_cpu = [0.057405134084524315
-    -0.16737701814667477
-     3.670092553395319
-     0.4052710265115234
-    -0.2342336895471848
-     0.07724866682808983
-     0.07379529136930803
-     0.004080775509334406
-    -0.12229701235980252
-    -0.0018175498192630844]
-
-    dual_infeas_4_cpu = [-0.0019510899995607155
-    0.07059765612753019
-    0.07104222896868044
-   -0.0061070561232362
-   -0.17080362240460334
-    0.0027843138173676155
-    0.003450443341928419
-    0.07157328724238488
-    0.07186014778306046
-   -0.28991023775276664]
-
-    dual_infeas_5_cpu = [-4.200422692029046
-    -0.33666635407314044
-     0.6407231013041463
-     0.0753794333813387
-     0.08217499064640721
-    -0.18387848500073856
-     1.4609510208883514
-     0.26889163927860954
-    -0.13184568446925882
-     0.065799538338272
-     0.04433557788772256
-     0.06271507036787678
-     2.807937471979002
-     0.0728632119983576
-    -0.20915786828282304
-     0.03875039932653859
-     0.020549534704044203]
-
-     dual_infeas_cpu = vcat(dual_infeas_1_cpu,dual_infeas_2_cpu,dual_infeas_3_cpu,dual_infeas_4_cpu,dual_infeas_5_cpu)
-
-     lambda_cpu = [0.488894     435.795        282.155     -0.0489912    248.787        348.203         0.135763     267.164        315.305
-     0.790211       0.88808        0.816501  -0.0653646      0.576725       0.456339      0.0687841      0.210784       0.140033
-    -8.67362e-19   -8.67362e-19   -0.0       -0.0           -4.33681e-19   -1.73472e-18  -0.0           -8.67362e-19   -8.67362e-19
-    -2.60209e-18   -0.0           -0.0       -1.73472e-18   -8.67362e-19   -8.67362e-19  -8.67362e-19   -0.0           -8.67362e-19]
-
-     @test norm(dpg_sol_cpu - Array(mod3.dpg_sol),  Inf) <= atol
-     @test norm(dqg_sol_cpu - Array(mod3.dqg_sol),  Inf) <= atol
-     @test norm(dline_var_cpu - Array(mod3.dline_var),  Inf) <= atol
-     @test norm(dline_fl_cpu - Array(mod3.dline_fl),  Inf) <= atol
-     @test norm(dtheta_sol_cpu - Array(mod3.dtheta_sol),  Inf) <= atol
-     @test norm(dw_sol_cpu - Array(mod3.dw_sol),  Inf) <= atol
-     @test norm(dual_infeas_cpu - Array(mod3.dual_infeas),  Inf) <= atol
-     @test norm((lambda_cpu - Array(mod3.lambda))./max.(abs.(lambda_cpu),1), Inf) <= 2e-6
+
+    dpg_sol_cpu = [-0.11550471288823375, 0.06549963131056268, 0.05367626157137195]
+    dqg_sol_cpu = [-0.005718588343597756, 0.01851317137973553, 0.01058363347077568]
+    dline_var_cpu = [-0.005569819650835843 -0.005054193193183796 -0.0031356731585941602 0.0003340424001857633 9.018007914709618e-6 -8.955231272232481e-6 0.0005809649468087098 -0.0035327541819175487 -0.0051673927876596205; -0.006653071419384663 -0.005109374028865239 -0.009433791106253716 0.003145428964346758 -0.0003118874694407943 -0.00022277448888749965 -0.004093726855025406 0.009746743401438716 0.005135828873903401; -0.0058992105835546605 -0.005767638775672478 -0.004596762982732082 0.0009542434497012447 -1.0647877468628053e-8 -1.492916208770773e-6 -1.068039123808262e-8 -1.068039123808262e-8 -0.004939127111986979; -0.005767628302723796 -0.004596740386993611 -1.0647877468628053e-8 -1.0647877468628053e-8 -1.4557418980203987e-6 -1.068039123808262e-8 0.0017380383383034281 -0.004939090417984399 -0.005767609350738688; -3.25491739361945e-19 0.005338385069326234 0.009542093395593967 0.0202615747235613 0.01764627076409719 0.01790665022472073 0.018092699885635936 0.018092713807391705 0.009562285691049546; 0.005338374345348201 0.009542068575890466 0.017646245194261424 0.017646257498616476 0.01790662375090709 0.018092686631961837 0.021440809620779835 0.009562257967646896 0.005338358797672126]
+    dline_fl_cpu = [-0.11550471214209485 -0.0550886536682177 -0.0545914661320966 0.05367626219021771 -0.003062020230877622 -0.0030394927033266204 -0.0654996296804065 0.06243360473855911 0.059908274544507235; -0.0057185925819239125 0.002880225309702243 0.004752099132814017 0.010583635657260775 0.00027191961424352743 0.0004625834867261366 -0.009295610035199166 0.009533344180015505 -0.003927910576173748; 0.11550471214209485 0.05459146673835657 0.056738283619246 -0.05367626219021771 0.003039493246316152 0.0030660246149609373 0.0654996296804065 -0.05990827415794294 -0.06041605732594946; -0.0034341779841658483 -0.004752098035100081 0.00542864572694711 -0.005700563960123412 -0.0004625831869150321 -0.00023773057022993016 0.018513174263915494 0.00392790777176107 0.0005539497426347278]
+    dtheta_sol_cpu = [-3.25491739361945e-19, 0.021440809620779835, 0.0202615747235613, 0.005338372737448853, 0.009542080985742217, 0.017646257818991696, 0.01790663698781391, 0.018092700108329824, 0.009562271829348221]
+    dw_sol_cpu = [-0.0058992105835546605, 0.0017380383383034281, 0.0009542434497012447, -0.005767625476378321, -0.004596751684862847, -1.0647877468628053e-8, -1.4743290533955858e-6, -1.068039123808262e-8, -0.004939108764985689]
+    dual_infeas_cpu = [-254.11036835411423, 111.34937322795656, 131.50684084986128, 0.3119989408890842, -8.034251796583263, -0.22792671306485607, -0.14821201491131938, 0.011898057499073345, 0.009455439882319458, 0.03537463279074429, -2.345388987728991, -0.18156150826834788, 0.08875478602088671, 0.03139112597637954, 0.02813068860448726, -0.11748530455855485, -1.2620894362848545, -0.10695004291180904, 0.19961093972437144, 0.05088261867417226, 0.057870735685249315, -0.1686020205956723, 3.6700753782985167, 0.3650331971871314, -0.1936042960743483, 0.0776146861999613, 0.0740166426887498, 0.004061310842689346, -0.1217101785537621, -0.0018110619226406265, -0.0019394705261958892, 0.0708845572940984, 0.07132702076591872, -0.006077871401322642, -0.169983866840654, 0.002767546200533025, 0.0034373043082680165, 0.0718559407428788, 0.07214140668385148, -0.29046051526063166, -4.20040217042528, -0.2957215255096348, 0.5993737191181183, 0.07574586733811475, 0.08238817068754835, -0.1824165829874273, 1.4673408402946906, 0.25436253640211093, -0.11631341488058167, 0.0658499676530561, 0.04476991944709829, 0.06513389823348271, 2.8151683733740307, 0.04757957575153994, -0.18020940646433298, 0.038415832071539036, 0.02118674588334765]
+    lambda_cpu = [0.4390310148492841 435.7732781204522 282.1544242381189 -0.020994217416346203 248.8021373471199 348.2250378203491 0.16468763831157662 267.1641217085404 315.2820765395896; 0.7925836307556712 0.8917093439553776 0.8196728616326903 -0.06555427114780897 0.5792170456138986 0.45857551396302043 0.06884954532846001 0.21271856587119875 0.14205944713641747; -1.734723475976807e-18 -8.673617379884035e-19 -2.168404344971009e-19 -0.0 -4.336808689942018e-19 -0.0 -8.673617379884035e-19 -8.673617379884035e-19 -8.673617379884035e-19; -1.734723475976807e-18 -0.0 -0.0 -0.0 -4.336808689942018e-19 -0.0 -1.734723475976807e-18 -0.0 -0.0]
+
+    @test norm(dpg_sol_cpu - Array(mod3.dpg_sol), Inf) <= atol
+    @test norm(dqg_sol_cpu - Array(mod3.dqg_sol), Inf) <= atol
+    @test norm(dline_var_cpu - Array(mod3.dline_var), Inf) <= atol
+    @test norm(dline_fl_cpu - Array(mod3.dline_fl), Inf) <= atol
+    @test norm(dtheta_sol_cpu - Array(mod3.dtheta_sol), Inf) <= atol
+    @test norm(dw_sol_cpu - Array(mod3.dw_sol), Inf) <= atol
+    @test norm(dual_infeas_cpu - Array(mod3.dual_infeas), Inf) <= atol
+    @test norm(lambda_cpu - Array(mod3.lambda), Inf) <= atol
 end
\ No newline at end of file