diff --git a/test/linear_program.jl b/test/linear_program.jl
new file mode 100644
index 00000000..aa90ed8d
--- /dev/null
+++ b/test/linear_program.jl
@@ -0,0 +1,250 @@
+# Copyright (c) 2020: Akshay Sharma and contributors
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+module TestLinearProgram
+
+using Test
+import DiffOpt
+import HiGHS
+import Ipopt
+import MathOptInterface as MOI
+import SCS
+
+const ATOL = 2e-4
+const RTOL = 2e-4
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$name", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+    return
+end
+
+include(joinpath(@__DIR__, "utils.jl"))
+
+function test_differentiating_LP_checking_gradients_for_non_active_contraints()
+    # Issue #40 from Gurobi.jl
+    # min  x
+    # s.t. x >= 0
+    #      x >= 3
+    nz = 1
+    qp_test_with_solutions(
+        HiGHS.Optimizer;
+        q = ones(nz),
+        G = -ones(2, nz),
+        h = [0.0, -3.0],
+        dzb = ones(nz),
+        dqf = ones(nz),
+        # Expected solutions
+        dGb = [0.0, 3.0],
+        dhb = [0.0, -1.0],
+    )
+    return
+end
+
+function test_differentiating_a_simple_LP_with_GreaterThan_constraint()
+    # this is canonically same as above test
+    # min  x
+    # s.t. x >= 3
+    nz = 1
+    qp_test_with_solutions(
+        Ipopt.Optimizer;
+        q = ones(nz),
+        G = -ones(1, nz),
+        h = [-3.0],
+        dzb = ones(nz),
+        dqf = ones(nz),
+        # Expected solutions
+        dGb = [3.0],
+        dhb = [-1.0],
+    )
+    return
+end
+
+function test_differentiating_lp()
+    # refered from - https://en.wikipedia.org/wiki/Simplex_algorithm#Example
+    # max 2x + 3y + 4z
+    # s.t. 3x+2y+z <= 10
+    #      2x+5y+3z <= 15
+    #      x,y,z >= 0
+    nz = 3
+    qp_test_with_solutions(
+        SCS.Optimizer;
+        q = [-2.0, -3.0, -4.0],
+        G = [
+            3.0 2.0 1.0
+            2.0 5.0 3.0
+            -1.0 0.0 0.0
+            0.0 -1.0 0.0
+            0.0 0.0 -1.0
+        ],
+        h = [10.0, 15.0, 0.0, 0.0, 0.0],
+        dzb = ones(nz),
+        dqf = ones(nz),
+        # Expected solutions
+        dqb = zeros(nz),
+        dGb = [
+            0.0 0.0 0.0
+            0.0 0.0 -5/3
+            0.0 0.0 5/3
+            0.0 0.0 -10/3
+            0.0 0.0 0.0
+        ],
+        dhb = [0.0, 1 / 3, -1 / 3, 2 / 3, 0.0],
+    )
+    return
+end
+
+function test_differentiating_LP_with_variable_bounds()
+    # max 2x + 3y + 4z
+    # s.t. 3x+2y+z <= 10
+    #      2x+5y+3z <= 15
+    #      x ≤ 3
+    #      0 ≤ y ≤ 2
+    #      z ≥ 2
+    #      x,y,z >= 0
+    # variant of previous test with same solution
+    nz = 3
+    qp_test_with_solutions(
+        HiGHS.Optimizer;
+        q = [-2.0, -3.0, -4.0],
+        G = [
+            3.0 2.0 1.0
+            2.0 5.0 3.0
+            -1.0 0.0 0.0
+            0.0 -1.0 0.0
+            0.0 0.0 -1.0
+            1.0 0.0 0.0
+            0.0 1.0 0.0
+            0.0 0.0 1.0
+        ],
+        h = [10.0, 15.0, 0.0, 0.0, 0.0, 3.0, 2.0, 6.0],
+        dzb = ones(nz),
+        dqf = ones(nz),
+        # Expected solutions
+        dqb = zeros(nz),
+        dGb = [
+            0.0 0.0 0.0
+            0.0 0.0 -5/3
+            0.0 0.0 5/3
+            0.0 0.0 -10/3
+            0.0 0.0 0.0
+            0.0 0.0 0.0
+            0.0 0.0 0.0
+            0.0 0.0 0.0
+        ],
+        dhb = [0.0, 1 / 3, -1 / 3, 2 / 3, 0.0, 0.0, 0.0, 0.0],
+    )
+    return
+end
+
+function test_differentiating_LP_with_variable_bounds_2()
+    nz = 3
+    qp_test_with_solutions(
+        HiGHS.Optimizer;
+        q = [-2.0, -3.0, -4.0],
+        G = [
+            3.0 2.0 1.0
+            2.0 5.0 3.0
+            0.0 -1.0 0.0
+            0.0 0.0 -1.0
+        ],
+        h = [10.0, 15.0, 0.0, 0.0],
+        fix_indices = [1],
+        fix_values = [0.0],
+        dzb = ones(nz),
+        dqf = ones(nz),
+        # Expected solutions
+        dqb = zeros(nz),
+        dGb = [
+            0.0 0.0 0.0
+            0.0 0.0 -5/3
+            0.0 0.0 -10/3
+            0.0 0.0 0.0
+        ],
+        dhb = [0.0, 1 / 3, 2 / 3, 0.0],
+        dAb = [0.0 0.0 -5 / 3],
+        dbb = [1 / 3],
+    )
+    return
+end
+
+"""
+max  2x + 3y + 4z
+s.t. 3x+2y+z <= 10
+     2x+5y+3z <= 15
+     x ≤ 3
+     0 ≤ y ≤ 2
+     6 ≥ z ≥ -1
+     x, y, z >= 0
+variant of previous test with same solution
+"""
+function test_differentiating_lp_with_saf_with_le_ge_constraints()
+    nz = 3
+    qp_test_with_solutions(
+        HiGHS.Optimizer;
+        q = [-2.0, -3.0, -4.0],
+        G = [
+            3.0 2.0 1.0
+            2.0 5.0 3.0
+            -1.0 0.0 0.0
+            0.0 -1.0 0.0
+            0.0 0.0 -1.0
+        ],
+        h = [10.0, 15.0, 0.0, 0.0, 0.0], #5
+        ub_indices = [1, 2, 3],
+        ub_values = [3.0, 2.0, 6.0],
+        lb_indices = [1],
+        lb_values = [-1.0],
+        dzb = ones(nz),
+        dqb = zeros(nz),
+        dGb = [
+            0.0 0.0 0.0
+            0.0 0.0 -5/3
+            0.0 0.0 5/3
+            0.0 0.0 -10/3
+            0.0 0.0 0.0
+            0.0 0.0 0.0
+            0.0 0.0 0.0
+            0.0 0.0 0.0
+            0.0 0.0 0.0
+        ],
+        dhb = [0.0, 1 / 3, -1 / 3, 2 / 3, 0.0, 0.0, 0.0, 0.0, 0.0],
+    )
+    return
+end
+
+function test_differentiating_lp_with_nonactive_constraints()
+    # Issue #40 from Gurobi.jl
+    # min  x
+    # s.t. x >= 0
+    #      x >= 3
+    qp_test_with_solutions(
+        HiGHS.Optimizer;
+        q = [1.0],
+        G = -ones(2, 1),
+        h = [0.0, -3.0],
+        dzb = -ones(1),
+        dhf = [0.0, 1.0],
+        # Expected solutions
+        z = [3.0],
+        λ = [0.0, 1.0],
+        dzf = -ones(1),
+        dλb = zeros(2),
+        dhb = [0.0, 1.0],
+        ∇zb = zeros(1),
+        ∇λb = [0.0, -1.0],
+        dλf = zeros(2),
+    )
+    return
+end
+
+end  # module
+
+TestLinearProgram.runtests()
diff --git a/test/quadratic_program.jl b/test/quadratic_program.jl
index 2c3f4048..4e535835 100644
--- a/test/quadratic_program.jl
+++ b/test/quadratic_program.jl
@@ -348,223 +348,6 @@ function test_differentiating_non_trivial_convex_qp_moi()
     return
 end
 
-function test_differentiating_LP_checking_gradients_for_non_active_contraints()
-    # Issue #40 from Gurobi.jl
-    # min  x
-    # s.t. x >= 0
-    #      x >= 3
-    nz = 1
-    qp_test_with_solutions(
-        HiGHS.Optimizer;
-        q = ones(nz),
-        G = -ones(2, nz),
-        h = [0.0, -3.0],
-        dzb = ones(nz),
-        dqf = ones(nz),
-        # Expected solutions
-        dGb = [0.0, 3.0],
-        dhb = [0.0, -1.0],
-    )
-    return
-end
-
-function test_differentiating_a_simple_LP_with_GreaterThan_constraint()
-    # this is canonically same as above test
-    # min  x
-    # s.t. x >= 3
-    nz = 1
-    qp_test_with_solutions(
-        Ipopt.Optimizer;
-        q = ones(nz),
-        G = -ones(1, nz),
-        h = [-3.0],
-        dzb = ones(nz),
-        dqf = ones(nz),
-        # Expected solutions
-        dGb = [3.0],
-        dhb = [-1.0],
-    )
-    return
-end
-
-function test_differentiating_lp()
-    # refered from - https://en.wikipedia.org/wiki/Simplex_algorithm#Example
-    # max 2x + 3y + 4z
-    # s.t. 3x+2y+z <= 10
-    #      2x+5y+3z <= 15
-    #      x,y,z >= 0
-    nz = 3
-    qp_test_with_solutions(
-        SCS.Optimizer;
-        q = [-2.0, -3.0, -4.0],
-        G = [
-            3.0 2.0 1.0
-            2.0 5.0 3.0
-            -1.0 0.0 0.0
-            0.0 -1.0 0.0
-            0.0 0.0 -1.0
-        ],
-        h = [10.0, 15.0, 0.0, 0.0, 0.0],
-        dzb = ones(nz),
-        dqf = ones(nz),
-        # Expected solutions
-        dqb = zeros(nz),
-        dGb = [
-            0.0 0.0 0.0
-            0.0 0.0 -5/3
-            0.0 0.0 5/3
-            0.0 0.0 -10/3
-            0.0 0.0 0.0
-        ],
-        dhb = [0.0, 1 / 3, -1 / 3, 2 / 3, 0.0],
-    )
-    return
-end
-
-function test_differentiating_LP_with_variable_bounds()
-    # max 2x + 3y + 4z
-    # s.t. 3x+2y+z <= 10
-    #      2x+5y+3z <= 15
-    #      x ≤ 3
-    #      0 ≤ y ≤ 2
-    #      z ≥ 2
-    #      x,y,z >= 0
-    # variant of previous test with same solution
-    nz = 3
-    qp_test_with_solutions(
-        HiGHS.Optimizer;
-        q = [-2.0, -3.0, -4.0],
-        G = [
-            3.0 2.0 1.0
-            2.0 5.0 3.0
-            -1.0 0.0 0.0
-            0.0 -1.0 0.0
-            0.0 0.0 -1.0
-            1.0 0.0 0.0
-            0.0 1.0 0.0
-            0.0 0.0 1.0
-        ],
-        h = [10.0, 15.0, 0.0, 0.0, 0.0, 3.0, 2.0, 6.0],
-        dzb = ones(nz),
-        dqf = ones(nz),
-        # Expected solutions
-        dqb = zeros(nz),
-        dGb = [
-            0.0 0.0 0.0
-            0.0 0.0 -5/3
-            0.0 0.0 5/3
-            0.0 0.0 -10/3
-            0.0 0.0 0.0
-            0.0 0.0 0.0
-            0.0 0.0 0.0
-            0.0 0.0 0.0
-        ],
-        dhb = [0.0, 1 / 3, -1 / 3, 2 / 3, 0.0, 0.0, 0.0, 0.0],
-    )
-    return
-end
-
-function test_differentiating_LP_with_variable_bounds_2()
-    nz = 3
-    qp_test_with_solutions(
-        HiGHS.Optimizer;
-        q = [-2.0, -3.0, -4.0],
-        G = [
-            3.0 2.0 1.0
-            2.0 5.0 3.0
-            0.0 -1.0 0.0
-            0.0 0.0 -1.0
-        ],
-        h = [10.0, 15.0, 0.0, 0.0],
-        fix_indices = [1],
-        fix_values = [0.0],
-        dzb = ones(nz),
-        dqf = ones(nz),
-        # Expected solutions
-        dqb = zeros(nz),
-        dGb = [
-            0.0 0.0 0.0
-            0.0 0.0 -5/3
-            0.0 0.0 -10/3
-            0.0 0.0 0.0
-        ],
-        dhb = [0.0, 1 / 3, 2 / 3, 0.0],
-        dAb = [0.0 0.0 -5 / 3],
-        dbb = [1 / 3],
-    )
-    return
-end
-
-"""
-max  2x + 3y + 4z
-s.t. 3x+2y+z <= 10
-     2x+5y+3z <= 15
-     x ≤ 3
-     0 ≤ y ≤ 2
-     6 ≥ z ≥ -1
-     x, y, z >= 0
-variant of previous test with same solution
-"""
-function test_differentiating_lp_with_saf_with_le_ge_constraints()
-    nz = 3
-    qp_test_with_solutions(
-        HiGHS.Optimizer;
-        q = [-2.0, -3.0, -4.0],
-        G = [
-            3.0 2.0 1.0
-            2.0 5.0 3.0
-            -1.0 0.0 0.0
-            0.0 -1.0 0.0
-            0.0 0.0 -1.0
-        ],
-        h = [10.0, 15.0, 0.0, 0.0, 0.0], #5
-        ub_indices = [1, 2, 3],
-        ub_values = [3.0, 2.0, 6.0],
-        lb_indices = [1],
-        lb_values = [-1.0],
-        dzb = ones(nz),
-        dqb = zeros(nz),
-        dGb = [
-            0.0 0.0 0.0
-            0.0 0.0 -5/3
-            0.0 0.0 5/3
-            0.0 0.0 -10/3
-            0.0 0.0 0.0
-            0.0 0.0 0.0
-            0.0 0.0 0.0
-            0.0 0.0 0.0
-            0.0 0.0 0.0
-        ],
-        dhb = [0.0, 1 / 3, -1 / 3, 2 / 3, 0.0, 0.0, 0.0, 0.0, 0.0],
-    )
-    return
-end
-
-function test_differentiating_lp_with_nonactive_constraints()
-    # Issue #40 from Gurobi.jl
-    # min  x
-    # s.t. x >= 0
-    #      x >= 3
-    qp_test_with_solutions(
-        HiGHS.Optimizer;
-        q = [1.0],
-        G = -ones(2, 1),
-        h = [0.0, -3.0],
-        dzb = -ones(1),
-        dhf = [0.0, 1.0],
-        # Expected solutions
-        z = [3.0],
-        λ = [0.0, 1.0],
-        dzf = -ones(1),
-        dλb = zeros(2),
-        dhb = [0.0, 1.0],
-        ∇zb = zeros(1),
-        ∇λb = [0.0, -1.0],
-        dλf = zeros(2),
-    )
-    return
-end
-
 end  # module
 
 TestQuadraticProgram.runtests()