diff --git a/netlib/getall.sh b/netlib/getall.sh
new file mode 100755
index 0000000..2acf63a
--- /dev/null
+++ b/netlib/getall.sh
@@ -0,0 +1,112 @@
+#!/bin/bash
+
+instances=(
+    "25fv47"
+    "80bau3b"
+    "adlittle"
+    "afiro"
+    "agg"
+    "agg2"
+    "agg3"
+    "ascii"
+    "bandm"
+    "beaconfd"
+    "blend"
+    "bnl1"
+    "bnl2"
+    "boeing1"
+    "boeing2"
+    "bore3d"
+    "brandy"
+    "capri"
+    "changes"
+    "cycle"
+    "czprob"
+    "d2q06c"
+    "d6cube"
+    "degen2"
+    "degen3"
+    "dfl001"
+    "e226"
+    "etamacro"
+    "fffff800"
+    "finnis"
+    "fit1d"
+    "fit1p"
+    "fit2d"
+    "fit2p"
+    "forplan"
+    "ganges"
+    "gfrd-pnc"
+    "greenbea"
+    "greenbeb"
+    "grow15"
+    "grow22"
+    "grow7"
+    "israel"
+    "kb2"
+    "kennington"
+    "lotfi"
+    "maros"
+    "maros-r7"
+    "minos"
+    "modszk1"
+    "nesm"
+    "perold"
+    "pilot"
+    "pilot.ja"
+    "pilot.we"
+    "pilot4"
+    "pilot87"
+    "pilotnov"
+    "recipe"
+    "sc105"
+    "sc205"
+    "sc50a"
+    "sc50b"
+    "scagr25"
+    "scagr7"
+    "scfxm1"
+    "scfxm2"
+    "scfxm3"
+    "scorpion"
+    "scrs8"
+    "scsd1"
+    "scsd6"
+    "scsd8"
+    "sctap1"
+    "sctap2"
+    "sctap3"
+    "seba"
+    "share1b"
+    "share2b"
+    "shell"
+    "ship04l"
+    "ship04s"
+    "ship08l"
+    "ship08s"
+    "ship12l"
+    "ship12s"
+    "sierra"
+    "stair"
+    "standata"
+    "standgub"
+    "standmps"
+    "stocfor1"
+    "stocfor2"
+    "stocfor3"
+    "stocfor3.old"
+    "truss"
+    "tuff"
+    "vtp.base"
+    "wood1p"
+    "woodw"
+)
+
+for i in "${instances[@]}"
+do
+    wget https://www.netlib.org/lp/data/$i
+    ./emps $i > $i.mps
+    rm $i
+done
+
diff --git a/src/Constant.jl b/src/Constant.jl
deleted file mode 100644
index 201db68..0000000
--- a/src/Constant.jl
+++ /dev/null
@@ -1,36 +0,0 @@
-"""
-Define constant variables
-"""
-
-const USE_INVB = true
-
-const INF = 1.0e+20
-const EPS = 1e-6
-
-const MAX_ITER = 30000000000
-const MAX_HISTORY = 100
-const PRINT_ITER_FREQ = 10
-
-@enum(
-    Status,
-    NotSolved,
-    Solve,
-    Optimal,
-    Unbounded,
-    Infeasible,
-    Feasible,
-    Cycle,
-)
-
-const BASIS_BASIC = 1
-const BASIS_AT_LOWER = 2
-const BASIS_AT_UPPER = 3
-const BASIS_FREE = 4
-
-@enum(
-    PivotRule,
-    Bland,
-    Steepest,
-    Dantzig,
-    Artificial,
-)
diff --git a/src/LP.jl b/src/LP.jl
index 26410f4..297910a 100644
--- a/src/LP.jl
+++ b/src/LP.jl
@@ -43,7 +43,7 @@ function canonical_form(standard::MatOI.LPForm{T, AT, VT}) where {T, AT, VT}
     @assert length(xl) == lp_ncols
     @assert length(xu) == lp_ncols
     for i = 1:nineq
-        if standard.c_lb[ineq[i]] > -INF
+        if isFinite(standard.c_lb[ineq[i]])
             xu[ncols(standard) + i] = standard.c_ub[ineq[i]] - standard.c_lb[ineq[i]]
         end
     end
@@ -54,10 +54,10 @@ function canonical_form(standard::MatOI.LPForm{T, AT, VT}) where {T, AT, VT}
     # create the submatrix for slack
     I = Int64[]; J = Int64[]; V = Float64[]; j = 1;
     for i in ineq
-        if standard.c_lb[i] > -INF
+        if isFinite(standard.c_lb[i])
             push!(I, i); push!(J, j); push!(V, -1.0);
             j += 1
-        elseif standard.c_ub[i] < INF
+        elseif isFinite(standard.c_ub[i])
             b[i] = standard.c_ub[i]
             push!(I, i); push!(J, j); push!(V, 1.0);
             j += 1
diff --git a/src/Parameters.jl b/src/Parameters.jl
new file mode 100644
index 0000000..c803289
--- /dev/null
+++ b/src/Parameters.jl
@@ -0,0 +1,10 @@
+"""
+    Parameters
+"""
+Base.@kwdef mutable struct Parameters
+    use_invB::Bool = true
+    max_iter::Int = 1000000000
+    max_history::Int = 100
+    print_iter_freq::Int = 10
+    pivot_rule::Type = DantzigRule
+end
diff --git a/src/PhaseOne/Artificial.jl b/src/PhaseOne/Artificial.jl
index 7de1d7d..8900cf2 100644
--- a/src/PhaseOne/Artificial.jl
+++ b/src/PhaseOne/Artificial.jl
@@ -143,7 +143,7 @@ function run(prob::MatOI.LPSolverForm, Tv::Type; kwargs...)::Simplex.SpxData
 
     # load the problem
     spx = Simplex.SpxData(p1lp, Tv)
-    spx.pivot_rule = deepcopy(pivot_rule)
+    spx.params.pivot_rule = pivot_rule
 
     # set basis
     basic = collect((ncols(p1lp)-nrows(p1lp)+1):ncols(p1lp))
diff --git a/src/PhaseOne/Cplex.jl b/src/PhaseOne/Cplex.jl
index f4ff3e0..41b9fb1 100644
--- a/src/PhaseOne/Cplex.jl
+++ b/src/PhaseOne/Cplex.jl
@@ -75,7 +75,7 @@ function run(lp::MatOI.LPSolverForm, Tv::Type; kwargs...)::Simplex.SpxData
 
     # load the problem
     spx = Simplex.SpxData(cpxlp, Tv)
-    spx.pivot_rule = deepcopy(pivot_rule)
+    spx.params.pivot_rule = deepcopy(pivot_rule)
 
     # set basis
     Simplex.set_basis(spx, B)
diff --git a/src/PivotRules.jl b/src/PivotRules.jl
index 2c8a0b7..fccb764 100644
--- a/src/PivotRules.jl
+++ b/src/PivotRules.jl
@@ -1,16 +1,17 @@
-# Bland's rule
-function pivot_Bland(spx::SpxData)
-    spx.enter = -1
+"""
+    BlandRule
 
-    # compute reduced cost
-    compute_reduced_cost(spx)
-    # @show spx.r
-    # @show spx.nonbasic
+Empty struct for Bland's pivot rule that can avoid cycling in pivoting
+"""
+struct BlandRule{T} <: AbstractPivotRule{T} end
 
+"Find an entering variable based on Bland's rule"
+function find_entering_variable!(::Type{BlandRule}, spx::SpxData)
+    spx.enter = -1
     min_enter = ncols(spx.lpdata) + 1
     for j = 1:length(spx.r)
         if spx.nonbasic[j] < min_enter
-            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
+            if spx.basis_status[spx.nonbasic[j]] == Basis_At_Upper
                 if spx.r[j] > 1e-10
                     spx.enter = spx.nonbasic[j]
                     spx.enter_pos = j
@@ -28,77 +29,91 @@ function pivot_Bland(spx::SpxData)
     # @show spx.enter, spx.r[spx.enter_pos]
 end
 
-# Steepest-edge rule
 """
-Steepest-edge rule:
+    SteepestEdge
+
+Mutable struct for steepest-edge pivoting rule:
 
 This is based on the following references.
 1. Bieling et al. An efficient GPU implementation of the revised simplex method. 2010
 2. Forest and Goldfarb. Steepest-edge simplex algorithms for linear programming. 1992
 3. Goldfarb and Reid. A practicable steepest-edge simplex algorithm. 1977
 """
-function pivot_steepest_edge(spx::SpxData)
-
-    compute_reduced_cost(spx)
-    # @show spx.r
-    # @show spx.nonbasic
 
-    compute_steepest_edge_weights(spx)
-
-    # compute the local slope
-    # @show length(spx.r), length(spx.gamma)
-    spx.s .= (spx.r).^2 ./ spx.gamma
-
-    spx.enter = -1
-    max_s = 0.0
-    for j = 1:length(spx.r)
-        if spx.s[j] > max_s
-            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER && sign(spx.r[j]) > 0
-                spx.enter = spx.nonbasic[j]
-                spx.enter_pos = j
-                max_s = spx.s[j]
-            elseif spx.basis_status[spx.nonbasic[j]] == BASIS_AT_LOWER && sign(spx.r[j]) < 0
-                spx.enter = spx.nonbasic[j]
-                spx.enter_pos = j
-                max_s = spx.s[j]
-            end
-        end
+mutable struct SteepestEdgeRule{T} <: AbstractPivotRule{T}
+    gamma::T # steepest-edge weights
+    v::T     # B^{-T} d
+    s::T     # = r.^2 ./ gamma
+    g::T     # invB * A_j
+    function SteepestEdgeRule(T::Type, m::Int, n::Int)
+        rule = new{T}()
+        rule.gamma = T{Float64}(undef, n - m)
+        rule.v = T{Float64}(undef, m)
+        rule.s = T{Float64}(undef, n - m)
+        rule.g = T{Float64}(undef, m)
+        return rule
     end
-    # @show spx.enter, s[spx.enter_pos]
 end
 
-function compute_steepest_edge_weights(spx::SpxData)
+"Find an entering variable based on steepest edge pivot rule"
+function find_entering_variable!(::Type{SteepestEdgeRule}, spx::SpxData)
+
+    pd = spx.pivot_data
     @timeit TO "compute steepest edge weights" begin
         if spx.iter == 1
             for j = 1:length(spx.nonbasic)
-                spx.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[j]]
-                spx.gamma[j] = spx.g' * spx.g
+                pd.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[j]]
+                pd.gamma[j] = pd.g' * pd.g
             end
         else
             # Here, spx.enter and spx.leave were obtained from the previous iteration.
             @timeit TO "gamma_j" begin
-                spx.gamma .+= 2 .* (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]) .* (spx.lpdata.A[:,spx.nonbasic]' * spx.v) .+ (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]).^2 .* (1 + spx.d' * spx.d)
+                pd.gamma .+= 2 .* (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]) .* (spx.lpdata.A[:,spx.nonbasic]' * pd.v) .+ (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]).^2 .* (1 + spx.d' * spx.d)
             end
 
             @timeit TO "gamm_e" begin
-                spx.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[spx.enter_pos]]
-                spx.gamma[spx.enter_pos] = spx.g' * spx.g
+                pd.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[spx.enter_pos]]
+                pd.gamma[spx.enter_pos] = pd.g' * pd.g
             end
         end
     end
-end
 
-# Dantzig's rule
-function pivot_Dantzig(spx::SpxData)
+    # compute the local slope
+    # @show length(spx.r), length(pd.gamma)
+    pd.s .= (spx.r).^2 ./ pd.gamma
+
     spx.enter = -1
+    max_s = 0.0
+    for j = 1:length(spx.r)
+        if pd.s[j] > max_s
+            if spx.basis_status[spx.nonbasic[j]] == Basis_At_Upper && sign(spx.r[j]) > 0
+                spx.enter = spx.nonbasic[j]
+                spx.enter_pos = j
+                max_s = pd.s[j]
+            elseif spx.basis_status[spx.nonbasic[j]] == Basis_At_Lower && sign(spx.r[j]) < 0
+                spx.enter = spx.nonbasic[j]
+                spx.enter_pos = j
+                max_s = pd.s[j]
+            end
+        end
+    end
+    # @show spx.enter, s[spx.enter_pos]
+end
+
+"""
+    DantzigRule
 
-    # compute reduced cost
-    compute_reduced_cost(spx)
+Struct for Dantzig's pivot rule
+"""
+struct DantzigRule{T} <: AbstractPivotRule{T} end
 
+"Find an entering variable based on Dantzig's pivot rule"
+function find_entering_variable!(::Type{DantzigRule}, spx::SpxData)
+    spx.enter = -1
     max_r = 0.0
     for j = 1:length(spx.r)
         if abs(spx.r[j]) > max_r
-            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
+            if spx.basis_status[spx.nonbasic[j]] == Basis_At_Upper
                 if spx.r[j] > 1e-6
                     spx.enter = spx.nonbasic[j]
                     spx.enter_pos = j
@@ -115,30 +130,3 @@ function pivot_Dantzig(spx::SpxData)
     end
     # @show spx.enter, spx.r[spx.enter_pos]
 end
-
-# Pivot to remove artificial variables from basis (based on Bland's rule)
-function pivot_to_remove_artificials(spx::SpxData)
-    spx.enter = -1
-
-    # compute reduced cost
-    compute_reduced_cost(spx)
-
-    min_enter = spx.start_artvars
-    for j = 1:length(spx.r)
-        if spx.nonbasic[j] < min_enter
-            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
-                if spx.r[j] > -1e-10
-                    spx.enter = spx.nonbasic[j]
-                    spx.enter_pos = j
-                    min_enter = spx.enter
-                end
-            else
-                if spx.r[j] < +1e-10
-                    spx.enter = spx.nonbasic[j]
-                    spx.enter_pos = j
-                    min_enter = spx.enter
-                end
-            end
-        end
-    end
-end
diff --git a/src/Simplex.jl b/src/Simplex.jl
index 7a683dc..10001a9 100644
--- a/src/Simplex.jl
+++ b/src/Simplex.jl
@@ -15,7 +15,9 @@ const TO = TimerOutput()
 nrows(lp::MatOI.AbstractLPForm{T}) where T = size(lp.A, 1)
 ncols(lp::MatOI.AbstractLPForm{T}) where T = size(lp.A, 2)
 
-include("Constant.jl")
+include("types.jl")
+include("Parameters.jl")
+include("Status.jl")
 include("GLPK.jl")
 include("LP.jl")
 include("PhaseOne/PhaseOne.jl")
@@ -23,7 +25,7 @@ include("PhaseOne/PhaseOne.jl")
 mutable struct SpxData{T<:AbstractArray}
     lpdata::MatOI.AbstractLPForm
 
-    basis_status::Vector{Int}
+    basis_status::Vector{BasisStatus}
     basic::Vector{Int}
     nonbasic::Vector{Int}
 
@@ -55,21 +57,15 @@ mutable struct SpxData{T<:AbstractArray}
     leave::Int
 
     update::Bool # whether to update basis or not
-    status::Status
+    status::AlgorithmStatus
 
     iter::Int
     history::Vector{Vector{Int}}
 
-    pivot_rule::PivotRule
+    params::Parameters
+    pivot_data::AbstractPivotRule{T}
 
-    # steepest-edge pivot only
-    gamma::T # steepest-edge weights
-    v::T     # B^{-T} d
-    s::T     # = r.^2 ./ gamma
-    g::T     # invB * A_j
-
-    # start index of the artificial variables
-    start_artvars::Int
+    array_type::Type
 
     function SpxData(lpdata::MatOI.AbstractLPForm{S}, T::Type) where S
         spx = new{T}()
@@ -79,6 +75,7 @@ mutable struct SpxData{T<:AbstractArray}
         if !in(T,[CuArray,Array])
             error("Unkown array type $(T).")
         end
+        spx.array_type = T
 
         # spx.view_basic_A = TArray{T,2}(undef, nrows(lpdata), nrows(lpdata))
         # spx.view_nonbasic_A = T{Float64,2}(undef, nrows(lpdata), ncols(lpdata) - nrows(lpdata))
@@ -103,11 +100,7 @@ mutable struct SpxData{T<:AbstractArray}
         spx.tl = T{Float64}(undef, nrows(lpdata))
         spx.tu = T{Float64}(undef, nrows(lpdata))
 
-        # Steepest-edge pivot only
-        spx.gamma = T{Float64}(undef, ncols(lpdata) - nrows(lpdata))
-        spx.v = T{Float64}(undef, nrows(lpdata))
-        spx.s = T{Float64}(undef, ncols(lpdata) - nrows(lpdata))
-        spx.g = T{Float64}(undef, nrows(lpdata))
+        spx.params = Parameters()
 
         fill!(spx.matQ, 0)
         fill!(spx.x, 0)
@@ -117,8 +110,6 @@ mutable struct SpxData{T<:AbstractArray}
         spx.iter = 1
         spx.history = Vector{Vector{Int}}()
 
-        spx.pivot_rule = Dantzig
-
         return spx
     end
 end
@@ -128,22 +119,14 @@ rhs(spx::SpxData) = spx.lpdata.b
 
 function inverse(spx::SpxData)
     @timeit TO "inverse" begin
-        # spx.invB .= spx.view_basic_A \ spx.matI
-        # @show size(spx.invB), size(spx.lpdata.A), size(spx.matI)
         spx.invB .= spx.lpdata.A[:,spx.basic] \ spx.matI
-        # droptol!(spx.invB, 1e-10)
     end
 end
 
 function update_inverse_basis(spx::SpxData)
     if spx.iter % 100 > 0
         @timeit TO "PF" begin
-            # view_invB_leave = @view spx.invB[spx.leave,:]
-            # view_matQ_leave = @view spx.matQ[spx.leave,:]
-
             # compute Q
-            # spx.matQ .= spx.d * view_invB_leave' ./ spx.d[spx.leave]
-            # view_matQ_leave .= view_invB_leave .* (1.0 - 1.0 / spx.d[spx.leave])
             spx.matQ .= spx.d * spx.invB[spx.leave,:]' ./ spx.d[spx.leave]
             spx.matQ[spx.leave,:] .= spx.invB[spx.leave,:] .* (1.0 - 1.0 / spx.d[spx.leave])
 
@@ -156,25 +139,23 @@ function update_inverse_basis(spx::SpxData)
     # @show mean(spx.invB), var(spx.invB), minimum(spx.invB), maximum(spx.invB)
 end
 
-function compute_xB(spx::SpxData)
-    @timeit TO "compute xB" begin
-        if USE_INVB
-            # spx.view_basic_x .= spx.invB * (rhs(spx) .- spx.view_nonbasic_A * spx.view_nonbasic_x)
-            spx.x[spx.basic] .= spx.invB * (rhs(spx) .- spx.lpdata.A[:,spx.nonbasic] * spx.x[spx.nonbasic])
-        else
-            # spx.view_basic_x .= spx.lpdata.A[:,spx.basic] \ (rhs(spx) .- spx.view_nonbasic_A * spx.view_nonbasic_x)
-            spx.x[spx.basic] .= spx.lpdata.A[:,spx.basic] \ (rhs(spx) .- spx.lpdata.A[:,spx.nonbasic] * spx.x[spx.nonbasic])
-        end
-        # @show spx.basic
-        # @show spx.x[spx.basic] .- spx.lpdata.v_lb[spx.basic]
-        @assert spx.x[spx.basic] >= spx.lpdata.v_lb[spx.basic]
-        @assert spx.x[spx.basic] <= spx.lpdata.v_ub[spx.basic]
+"Compute basic solution"
+function compute_xB!(spx::SpxData)
+    if spx.params.use_invB
+        # spx.view_basic_x .= spx.invB * (rhs(spx) .- spx.view_nonbasic_A * spx.view_nonbasic_x)
+        spx.x[spx.basic] .= spx.invB * (rhs(spx) .- spx.lpdata.A[:,spx.nonbasic] * spx.x[spx.nonbasic])
+    else
+        # spx.view_basic_x .= spx.lpdata.A[:,spx.basic] \ (rhs(spx) .- spx.view_nonbasic_A * spx.view_nonbasic_x)
+        spx.x[spx.basic] .= spx.lpdata.A[:,spx.basic] \ (rhs(spx) .- spx.lpdata.A[:,spx.nonbasic] * spx.x[spx.nonbasic])
     end
+    @assert spx.x[spx.basic] >= spx.lpdata.v_lb[spx.basic]
+    @assert spx.x[spx.basic] <= spx.lpdata.v_ub[spx.basic]
 end
 
-function compute_pB(spx::SpxData)
-    @timeit TO "compute pB" begin
-        if USE_INVB
+"Compute dual multiplier"
+function compute_pB!(spx::SpxData)
+    if spx.update
+        if spx.params.use_invB
             # spx.pB .= transpose(spx.invB) * spx.view_basic_c
             spx.pB .= transpose(spx.invB) * spx.lpdata.c[spx.basic]
         else
@@ -183,113 +164,73 @@ function compute_pB(spx::SpxData)
     end
 end
 
-function compute_reduced_cost(spx::SpxData)
-    @timeit TO "compute reduced cost" begin
-        # spx.r .= spx.view_nonbasic_c .- spx.view_nonbasic_A' * spx.pB
-        spx.r .= spx.lpdata.c[spx.nonbasic] .- spx.lpdata.A[:,spx.nonbasic]' * spx.pB
-    end
-end
-
-function compute_direction(spx::SpxData)
-    @timeit TO "compute direction" begin
-        if USE_INVB
-            spx.d .= spx.invB * spx.lpdata.A[:,spx.enter]
-        else
-            view_A_enter = @view spx.lpdata.A[:,spx.enter]
-            # @show cond(Matrix(spx.lpdata.A[:,spx.basic]))
-            spx.d .= spx.lpdata.A[:,spx.basic] \ view_A_enter
-        end
-    end
+"Compute reduced cost"
+function compute_reduced_cost!(spx::SpxData)
+    # spx.r .= spx.view_nonbasic_c .- spx.view_nonbasic_A' * spx.pB
+    spx.r .= spx.lpdata.c[spx.nonbasic] .- spx.lpdata.A[:,spx.nonbasic]' * spx.pB
 end
 
-function compute_entering_variable(spx::SpxData)
-    @timeit TO "compute entering variable" begin
-        if spx.pivot_rule == Bland
-            pivot_Bland(spx)
-        elseif spx.pivot_rule == Steepest
-            pivot_steepest_edge(spx)
-        elseif spx.pivot_rule == Dantzig
-            pivot_Dantzig(spx)
-        elseif spx.pivot_rule == Artificial
-            pivot_to_remove_artificials(spx)
-        else
-            @error("Invalid pivot rule.")
-        end
+"Compute direction"
+function compute_direction!(spx::SpxData)
+    if spx.params.use_invB
+        spx.d .= spx.invB * spx.lpdata.A[:,spx.enter]
+    else
+        view_A_enter = @view spx.lpdata.A[:,spx.enter]
+        # @show cond(Matrix(spx.lpdata.A[:,spx.basic]))
+        spx.d .= spx.lpdata.A[:,spx.basic] \ view_A_enter
     end
 end
 
-function pivot(spx::SpxData)
+"perform a ratio test and update solution"
+function ratio_test!(spx::SpxData)
     spx.leave = -1
     spx.update = true
     best_t = 0.0
 
-    # compute the direction
-    compute_direction(spx)
-    # @show norm(spx.d)
-
-    # Terminate with unboundedness
-    if norm(spx.d) < EPS
-        return
-    end
-
-    # ratio test
-    @timeit TO "ratio test" begin
-        spx.tl .= (spx.x[spx.basic] .- spx.lpdata.v_lb[spx.basic]) ./ spx.d
-        spx.tu .= (spx.x[spx.basic] .- spx.lpdata.v_ub[spx.basic]) ./ spx.d
-        if spx.r[spx.enter_pos] < 0
-            best_t = Inf
-            for i = 1:nrows(spx.lpdata)
-                if spx.d[i] < -EPS && best_t > spx.tu[i]
-                    best_t = spx.tu[i]
-                    spx.leave = i
-                end
-                if spx.d[i] > +EPS && best_t > spx.tl[i]
-                    best_t = spx.tl[i]
-                    spx.leave = i
-                end
+    spx.tl .= (spx.x[spx.basic] .- spx.lpdata.v_lb[spx.basic]) ./ spx.d
+    spx.tu .= (spx.x[spx.basic] .- spx.lpdata.v_ub[spx.basic]) ./ spx.d
+    if spx.r[spx.enter_pos] < 0
+        best_t = Inf
+        for i = 1:nrows(spx.lpdata)
+            if isNegative(spx.d[i]) && best_t > spx.tu[i]
+                best_t = spx.tu[i]
+                spx.leave = i
             end
-            if best_t > spx.lpdata.v_ub[spx.enter] - spx.x[spx.enter] && best_t < Inf
-                best_t = spx.lpdata.v_ub[spx.enter] - spx.lpdata.v_lb[spx.enter]
-                spx.x[spx.enter] = spx.lpdata.v_ub[spx.enter]
-                spx.basis_status[spx.enter] = BASIS_AT_UPPER
-                spx.update = false
+            if isPositive(spx.d[i]) && best_t > spx.tl[i]
+                best_t = spx.tl[i]
+                spx.leave = i
             end
-        elseif spx.r[spx.enter_pos] > 0
-            best_t = -Inf
-            for i = 1:nrows(spx.lpdata)
-                if spx.d[i] < -EPS && best_t < spx.tl[i]
-                    best_t = spx.tl[i]
-                    spx.leave = i
-                end
-                if spx.d[i] > +EPS && best_t < spx.tu[i]
-                    best_t = spx.tu[i]
-                    spx.leave = i
-                end
-                # if spx.pivot_rule == Artificial && spx.basic[i] < spx.start_artvars
-                #     spx.tu[i] = -Inf
-                # end
+        end
+        if best_t > spx.lpdata.v_ub[spx.enter] - spx.x[spx.enter] && best_t < Inf
+            best_t = spx.lpdata.v_ub[spx.enter] - spx.lpdata.v_lb[spx.enter]
+            spx.x[spx.enter] = spx.lpdata.v_ub[spx.enter]
+            spx.basis_status[spx.enter] = Basis_At_Upper
+            spx.update = false
+        end
+    elseif spx.r[spx.enter_pos] > 0
+        best_t = -Inf
+        for i = 1:nrows(spx.lpdata)
+            if isNegative(spx.d[i]) && best_t < spx.tl[i]
+                best_t = spx.tl[i]
+                spx.leave = i
             end
-            if best_t < spx.lpdata.v_lb[spx.enter] - spx.x[spx.enter] && best_t > -Inf
-                best_t = spx.lpdata.v_lb[spx.enter] - spx.lpdata.v_ub[spx.enter]
-                spx.x[spx.enter] = spx.lpdata.v_lb[spx.enter]
-                spx.basis_status[spx.enter] = BASIS_AT_LOWER
-                spx.update = false
+            if isPositive(spx.d[i]) && best_t < spx.tu[i]
+                best_t = spx.tu[i]
+                spx.leave = i
             end
-        else
-            @error("The reduced cost is zero.")
         end
-        @assert spx.leave > 0
-
-        best_t = abs(best_t*spx.d[spx.leave]) < EPS ? 0.0 : best_t
-        # @show abs(best_t*spx.d[spx.leave])
-        # if abs(best_t) < Inf
-        #     best_t = abs(best_t*spx.d[spx.leave]) < EPS ? 0.0 : best_t
-        # else
-        #     spx.leave = -1
-        #     return
-        # end
-        # @assert abs(spx.d[spx.leave]) >= 1e-8
+        if best_t < spx.lpdata.v_lb[spx.enter] - spx.x[spx.enter] && best_t > -Inf
+            best_t = spx.lpdata.v_lb[spx.enter] - spx.lpdata.v_ub[spx.enter]
+            spx.x[spx.enter] = spx.lpdata.v_lb[spx.enter]
+            spx.basis_status[spx.enter] = Basis_At_Lower
+            spx.update = false
+        end
+    else
+        @error("The reduced cost is zero.")
     end
+    @assert spx.leave > 0
+
+    best_t = isZero(best_t*spx.d[spx.leave]) ? 0.0 : best_t
 
     # update solution
     # @show spx.x[spx.basic[spx.leave]], spx.lpdata.v_lb[spx.basic[spx.leave]], spx.lpdata.v_ub[spx.basic[spx.leave]]
@@ -306,95 +247,90 @@ function pivot(spx::SpxData)
     end
 end
 
-function detect_cycle(spx::SpxData)
-    @timeit TO "detect cycly" begin
-        if in(spx.basic, spx.history)
-            @info "Cycle is detected."
-            spx.status = Cycle
-        end
-        if length(spx.history) == MAX_HISTORY
-            pop!(spx.history)
-        end
-        push!(spx.history, deepcopy(spx.basic))
+"Detect pivot cycle"
+function detect_cycle!(spx::SpxData)
+    if in(spx.basic, spx.history)
+        @info "Cycle is detected."
+        spx.status = Cycle
+    end
+    if length(spx.history) == spx.params.max_history
+        pop!(spx.history)
     end
+    push!(spx.history, deepcopy(spx.basic))
 end
 
-function update_basis(spx::SpxData)
-    @assert spx.d[spx.leave] != 0
+function update_basis!(spx::SpxData)
 
-    # update basis status
-    @timeit TO "update basis" begin
-        if isapprox(spx.x[spx.basic[spx.leave]], spx.lpdata.v_lb[spx.basic[spx.leave]], atol=1e-6)
-            spx.basis_status[spx.basic[spx.leave]] = BASIS_AT_LOWER
-            spx.x[spx.basic[spx.leave]] = spx.lpdata.v_lb[spx.basic[spx.leave]]
-        elseif isapprox(spx.x[spx.basic[spx.leave]], spx.lpdata.v_ub[spx.basic[spx.leave]], atol=1e-6)
-            spx.basis_status[spx.basic[spx.leave]] = BASIS_AT_UPPER
-            spx.x[spx.basic[spx.leave]] = spx.lpdata.v_ub[spx.basic[spx.leave]]
-        else
-            spx.basis_status[spx.basic[spx.leave]] = BASIS_FREE
-        end
-        spx.basis_status[spx.enter] = BASIS_BASIC
+    if spx.params.pivot_rule == SteepestEdgeRule
+        # compute v before changing basis
+        spx.pivot_data.v .= transpose(spx.invB) * spx.d
+    end
 
-        # update basic/nonbasic indices
-        spx.nonbasic[spx.enter_pos] = deepcopy(spx.basic[spx.leave])
-        spx.basic[spx.leave] = deepcopy(spx.enter)
-        # @show spx.basic
+    # update basis status
+    if isapprox(spx.x[spx.basic[spx.leave]], spx.lpdata.v_lb[spx.basic[spx.leave]], atol=1e-6)
+        spx.basis_status[spx.basic[spx.leave]] = Basis_At_Lower
+        spx.x[spx.basic[spx.leave]] = spx.lpdata.v_lb[spx.basic[spx.leave]]
+    elseif isapprox(spx.x[spx.basic[spx.leave]], spx.lpdata.v_ub[spx.basic[spx.leave]], atol=1e-6)
+        spx.basis_status[spx.basic[spx.leave]] = Basis_At_Upper
+        spx.x[spx.basic[spx.leave]] = spx.lpdata.v_ub[spx.basic[spx.leave]]
+    else
+        spx.basis_status[spx.basic[spx.leave]] = Basis_Free
     end
+    spx.basis_status[spx.enter] = Basis_Basic
 
-    detect_cycle(spx)
+    # update basic/nonbasic indices
+    spx.nonbasic[spx.enter_pos] = deepcopy(spx.basic[spx.leave])
+    spx.basic[spx.leave] = deepcopy(spx.enter)
+    # @show spx.basic
 
-    if USE_INVB
+    if spx.params.use_invB
         update_inverse_basis(spx)
     end
 end
 
+"Take one interation"
 function iterate(spx::SpxData)
-    @timeit TO "One iteration" begin
-        if spx.iter % PRINT_ITER_FREQ == 0 && spx.pivot_rule != Artificial
-            println("Iteration $(spx.iter): objective $(objective(spx))")
-        end
-
-        if spx.update
-            compute_pB(spx)
-            # @show spx.pB
-        end
-
-        compute_entering_variable(spx)
+    if spx.iter % spx.params.print_iter_freq == 0
+        println("Iteration $(spx.iter): objective $(objective(spx))")
+    end
 
-        if spx.enter < 0
-            spx.status = Optimal
-            return
-        end
+    @timeit TO "compute dual multiplier" compute_pB!(spx)
+    @timeit TO "compute reduced cost" compute_reduced_cost!(spx)
+    @timeit TO "find entering variable" find_entering_variable!(spx.params.pivot_rule, spx)
 
-        pivot(spx)
+    # Optimal if no entering variable is found
+    if spx.enter < 0
+        spx.status = Optimal
+        return
+    end
 
-        if spx.leave < 0
-            # spx.status = spx.pivot_rule == Artificial ? Feasible : Unbounded
-            spx.status = Unbounded
-            return
-        end
+    @timeit TO "compute direction" compute_direction!(spx)
 
-        # println("Update basis? ", spx.update)
+    # Terminate with unboundedness
+    if isNonPositive(norm(spx.d))
+        spx.status = Unbounded
+        return
+    end
 
-        if spx.update
-            if spx.pivot_rule == Steepest
-                @timeit TO "compute v" begin
-                    # compute v before changing basis
-                    spx.v .= transpose(spx.invB) * spx.d
-                end
-            end
+    @timeit TO "ratio test" ratio_test!(spx)
 
-            # println("  Entering/Leaving variables: $(spx.enter), $(spx.basic[spx.leave])")
-            update_basis(spx)
-        end
+    # Unbounded if no leaving variable is found
+    if spx.leave < 0
+        spx.status = Unbounded
+        return
+    end
 
-        spx.iter += 1
+    if spx.update
+        @timeit TO "update basis" update_basis!(spx)
+        @timeit TO "detect pivot cycle" detect_cycle!(spx)
     end
+
+    spx.iter += 1
 end
 
 function run(prob::MatOI.LPForm{T, AT, VT};
     basis::Vector{Int} = Int[],
-    pivot_rule::PivotRule = Dantzig,
+    pivot_rule::Type = DantzigRule,
     phase_one_method::PhaseOne.Method = PhaseOne.CPLEX,
     gpu = false,
     performance_io::IO = stdout) where {T, AT, VT}
@@ -408,9 +344,7 @@ function run(prob::MatOI.LPForm{T, AT, VT};
             presolve_prob = presolve(prob)
         end
 
-        @timeit TO "scaling" begin
-            scaling!(presolve_prob)
-        end
+        @timeit TO "scaling" scaling!(presolve_prob)
 
         # convert the problem into a canonical form
         canonical = Simplex.canonical_form(presolve_prob)
@@ -451,13 +385,11 @@ function run(prob::MatOI.LPForm{T, AT, VT};
             else
                 println("Basis information is provided.")
                 spx = SpxData(processed_prob, TT)
-                spx.pivot_rule = pivot_rule
+                spx.params.pivot_rule = pivot_rule
                 set_basis(spx, basis)
             end
 
-            @timeit TO "Phase 2" begin
-                run_core(spx)
-            end
+            @timeit TO "Phase 2" run_core(spx)
 
             status = spx.status
             println("Phase 2 is done: status $(Symbol(spx.status))")
@@ -472,29 +404,35 @@ function run(prob::MatOI.LPForm{T, AT, VT};
 end
 
 function run_core(spx::SpxData)
-    if USE_INVB
+
+    # Initialize pivot data for steepest edge rule
+    if spx.params.pivot_rule == SteepestEdgeRule && !isdefined(spx, :pivot_data)
+        spx.pivot_data = SteepestEdgeRule(spx.array_type, nrows(spx.lpdata), ncols(spx.lpdata))
+    end
+
+    if spx.params.use_invB
         inverse(spx)
     end
     # @show spx.invB
 
-    compute_xB(spx)
+    @timeit TO "compute basic solution" compute_xB!(spx)
     # @show spx.x
 
     # main iterations
-    my_pivot_rule = spx.pivot_rule
-    while spx.status in [Solve, Cycle] && spx.iter < MAX_ITER
+    my_pivot_rule = spx.params.pivot_rule
+    while spx.status in [Solve, Cycle] && spx.iter < spx.params.max_iter
         # Use Bland's rule if cycle is detected
         if spx.status == Cycle
-            spx.pivot_rule = Bland
+            spx.params.pivot_rule = BlandRule
         end
-        iterate(spx)
-        spx.pivot_rule = my_pivot_rule
+        @timeit TO "One iteration" iterate(spx)
+        spx.params.pivot_rule = my_pivot_rule
     end
 end
 
 include("Presolve.jl")
 include("PivotRules.jl")
 include("WarmStart.jl")
-# include("utils.jl")
+include("utils.jl")
 
 end # module
diff --git a/src/Status.jl b/src/Status.jl
new file mode 100644
index 0000000..b041f6f
--- /dev/null
+++ b/src/Status.jl
@@ -0,0 +1,31 @@
+"""
+    AlgorithmStatus
+
+This indicates the status of algorithm.
+"""
+@enum(AlgorithmStatus,
+    NotSolved,
+    Solve,
+    Optimal,
+    Unbounded,
+    Infeasible,
+    Feasible,
+    Cycle,
+)
+
+"""
+    BasisStatus
+
+- `Basis_Basic`: basic variable
+- `Basis_At_Lower`: nonbasic variable at lower bound
+- `Basis_At_Upper`: nonbasic variable at upper bound
+- `Basis_Fixed`: nonbasic variable at the fixed bound
+- `Basis_Free`: nonbasic free variable
+"""
+@enum(BasisStatus,
+    Basis_Basic,
+    Basis_At_Lower,
+    Basis_At_Upper,
+    Basis_Fixed,
+    Basis_Free
+)
diff --git a/src/WarmStart.jl b/src/WarmStart.jl
index 5ea65e7..fcbabb0 100644
--- a/src/WarmStart.jl
+++ b/src/WarmStart.jl
@@ -2,7 +2,7 @@
 
 function set_basis(spx::SpxData, basis::Vector{Int})
     # Set basic and nonbasic indices
-    spx.basis_status = Vector{Int}(undef, ncols(spx.lpdata))
+    spx.basis_status = Vector{BasisStatus}(undef, ncols(spx.lpdata))
     spx.basic = basis
     spx.nonbasic = Int[]
     sizehint!(spx.nonbasic, ncols(spx.lpdata) - nrows(spx.lpdata))
@@ -12,20 +12,20 @@ function set_basis(spx::SpxData, basis::Vector{Int})
             axl = abs(spx.lpdata.v_lb[j])
             axu = abs(spx.lpdata.v_ub[j])
             if spx.lpdata.v_lb[j] > -Inf && axl <= axu
-                spx.basis_status[j] = BASIS_AT_LOWER
+                spx.basis_status[j] = Basis_At_Lower
             elseif spx.lpdata.v_ub[j] < Inf && axl > axu
-                spx.basis_status[j] = BASIS_AT_UPPER
+                spx.basis_status[j] = Basis_At_Upper
             else
-                spx.basis_status[j] = BASIS_FREE
+                spx.basis_status[j] = Basis_Free
             end
         else
-            spx.basis_status[j] = BASIS_BASIC
+            spx.basis_status[j] = Basis_Basic
         end
     end
     set_nonbasic(spx)
 end
 
-function set_basis_status(spx::SpxData, basis_status::Vector{Int})
+function set_basis_status(spx::SpxData, basis_status::Vector{BasisStatus})
     # Set basic and nonbasic indices
     spx.basis_status = basis_status
     spx.basic = Int[]
@@ -33,7 +33,7 @@ function set_basis_status(spx::SpxData, basis_status::Vector{Int})
     sizehint!(spx.basic, nrows(spx.lpdata))
     sizehint!(spx.nonbasic, ncols(spx.lpdata) - nrows(spx.lpdata))
     for j = 1:ncols(spx.lpdata)
-        if basis_status[j] == BASIS_BASIC
+        if basis_status[j] == Basis_Basic
             push!(spx.basic, j)
         else
             push!(spx.nonbasic, j)
@@ -42,7 +42,7 @@ function set_basis_status(spx::SpxData, basis_status::Vector{Int})
     set_nonbasic(spx)
 end
 
-function set_basis(spx::SpxData, basis::Vector{Int}, basis_status::Vector{Int})
+function set_basis(spx::SpxData, basis::Vector{Int}, basis_status::Vector{BasisStatus})
     # Set basic and nonbasic indices
     spx.basis_status = basis_status
     spx.basic = basis
@@ -57,23 +57,23 @@ function set_basis(spx::SpxData, basis::Vector{Int}, basis_status::Vector{Int})
 end
 
 function set_basis(spx::SpxData{T}, x::T) where T
-    spx.basis_status = Vector{Int}(undef, ncols(spx.lpdata))
+    spx.basis_status = Vector{BasisStatus}(undef, ncols(spx.lpdata))
     spx.basic = Int[]
     spx.nonbasic = Int[]
     sizehint!(spx.basic, nrows(spx.lpdata))
     sizehint!(spx.nonbasic, ncols(spx.lpdata) - nrows(spx.lpdata))
     for j = 1:ncols(spx.lpdata)
         if x[j] == spx.lpdata.v_lb[j]
-            spx.basis_status[j] = BASIS_AT_LOWER
+            spx.basis_status[j] = Basis_At_Lower
             push!(spx.nonbasic, j)
         elseif x[j] == spx.lpdata.v_ub[j]
-            spx.basis_status[j] = BASIS_AT_UPPER
+            spx.basis_status[j] = Basis_At_Upper
             push!(spx.nonbasic, j)
         elseif length(spx.basic) < nrows(spx.lpdata)
-            spx.basis_status[j] = BASIS_BASIC
+            spx.basis_status[j] = Basis_Basic
             push!(spx.basic, j)
         else
-            spx.basis_status[j] = BASIS_FREE
+            spx.basis_status[j] = Basis_Free
             push!(spx.nonbasic, j)
         end
     end
@@ -82,11 +82,11 @@ end
 
 function set_nonbasic(spx::SpxData)
     for j in spx.nonbasic
-        if spx.basis_status[j] == BASIS_AT_LOWER
+        if spx.basis_status[j] == Basis_At_Lower
             spx.x[j] = spx.lpdata.v_lb[j]
-        elseif spx.basis_status[j] == BASIS_AT_UPPER
+        elseif spx.basis_status[j] == Basis_At_Upper
             spx.x[j] = spx.lpdata.v_ub[j]
-        elseif spx.basis_status[j] == BASIS_FREE
+        elseif spx.basis_status[j] == Basis_Free
             # @warn("Free nonbasic variable (x$j) is set to zero.")
             spx.x[j] = 0.0
         end
diff --git a/src/types.jl b/src/types.jl
new file mode 100644
index 0000000..6358063
--- /dev/null
+++ b/src/types.jl
@@ -0,0 +1,6 @@
+"""
+    AbstractPivotRule
+
+Abstract type of pivot rules
+"""
+abstract type AbstractPivotRule{T} end
diff --git a/src/utils.jl b/src/utils.jl
index 0e7556e..41709b1 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,79 +1,16 @@
-function print(lp::StandardLpData)
-    println("Number of rows   : $(lp.nrows)")
-    println("Number of columns: $(lp.ncols)")
-    println("Constraint matrix: $(Matrix(lp.A))")
-    println("Row lower bounds: $(lp.bl)")
-    println("Row upper bounds: $(lp.bu)")
-    println("Objective coefficients: $(lp.c)")
-    println("Column lower bound: $(lp.xl)")
-    println("Column upper bound: $(lp.xu)")
-end
 
-function print(lp::CanonicalLpData)
-    println("Number of rows   : $(lp.nrows)")
-    println("Number of columns: $(lp.ncols)")
-    println("Constraint matrix: $(Matrix(lp.A))")
-    println("Row bounds: $(lp.b)")
-    println("Objective coefficients: $(lp.c)")
-    println("Column lower bound: $(lp.xl)")
-    println("Column upper bound: $(lp.xu)")
-end
+const INF = 1.0e+20
+const EPS = 1.0e-6
 
-function summary(lp::StandardLpData)
-    println("Number of rows   : $(lp.nrows)")
-    println("Number of columns: $(lp.ncols)")
-    # println("A: Min $(minimum(nonzeros(lp.A))), Max $(maximum(nonzeros(lp.A)))")
-    println("bl: Min $(minimum(lp.bl)), Max $(maximum(lp.bl))")
-    println("bu: Min $(minimum(lp.bu)), Max $(maximum(lp.bu))")
-    println("xl: Min $(minimum(lp.xl)), Max $(maximum(lp.xl))")
-    println("xu: Min $(minimum(lp.xu)), Max $(maximum(lp.xu))")
-    println("c: Min $(minimum(lp.c)), Max $(maximum(lp.c))")
-
-    cons = Dict("EQ"=>0, "GT"=>0, "LT"=>0, "RNG"=>0, "FR"=>0)
-    vars = Dict("EQ"=>0, "GT"=>0, "LT"=>0, "RNG"=>0, "FR"=>0)
-    for i = 1:lp.nrows
-        if lp.bl[i] == lp.bu[i]
-            cons["EQ"] += 1
-        elseif lp.bl[i] > -Inf
-            if lp.bu[i] < Inf
-                cons["RNG"] += 1
-            else
-                cons["GT"] += 1
-            end
-        elseif lp.bu[i] <= Inf
-            cons["LT"] += 1
-        else
-            cons["FR"] += 1
-        end
-    end
-
-    for j = 1:lp.ncols
-        if lp.xl[j] == lp.xu[j]
-            vars["EQ"] += 1
-        elseif lp.xl[j] > -Inf
-            if lp.xu[j] < Inf
-                vars["RNG"] += 1
-            else
-                vars["GT"] += 1
-            end
-        elseif lp.xu[j] <= Inf
-            vars["LT"] += 1
-        else
-            vars["FR"] += 1
-        end
-    end
-
-    println("Constraints:")
-    for (k,v) = cons
-        println(" $k = $v")
-    end
-    println("Variables:")
-    for (k,v) = vars
-        println(" $k = $v")
-    end
-
-    # for j = 1:lp.ncols
-    #     if lp.xl[j] 
-    # end
-    # println("Free: ")
-end
\ No newline at end of file
+isGT(x::Real, y::Real, tol::Real = EPS)::Bool = (x - y > +tol)
+isGE(x::Real, y::Real, tol::Real = EPS)::Bool = (x - y >= -tol)
+isLE(x::Real, y::Real, tol::Real = EPS)::Bool = (x - y <= +tol)
+isLT(x::Real, y::Real, tol::Real = EPS)::Bool = (x - y < -tol)
+isEQ(x::Real, y::Real, tol::Real = EPS)::Bool = (abs(x-y) <= tol)
+isPositive(x::Real, tol::Real = EPS)::Bool = (x > tol)
+isNonNegative(x::Real, tol::Real = EPS)::Bool = (x >= -tol)
+isNonPositive(x::Real, tol::Real = EPS)::Bool = (x < tol)
+isNegative(x::Real, tol::Real = EPS)::Bool = (x < -tol)
+isInf(x::Real, tol::Real = INF)::Bool = (x > tol)
+isFinite(x::Real, tol::Real = INF)::Bool = (abs(x) < tol)
+isZero(x::Real, tol::Real = EPS)::Bool = (abs(x) < EPS)
diff --git a/test/netlib.jl b/test/netlib.jl
index 03f282b..6519980 100644
--- a/test/netlib.jl
+++ b/test/netlib.jl
@@ -16,12 +16,12 @@ phaseone_methods = [
     Simplex.PhaseOne.CPLEX,
 ]
 pivot_rules = [
-    Simplex.Bland,
-    Simplex.Steepest,
-    Simplex.Dantzig,
+    Simplex.BlandRule,
+    Simplex.SteepestEdgeRule,
+    Simplex.DantzigRule,
 ]
 
-function run_netlib_instance(netlib, use_gpu, method, pivot)::Simplex.Status
+function run_netlib_instance(netlib, use_gpu, method, pivot)::Simplex.AlgorithmStatus
     c, xl, xu, bl, bu, A = Simplex.GLPK_read_mps(netlib)
     lp =  MatOI.LPForm{Float64, typeof(A), typeof(c)}(
         MOI.MIN_SENSE, c, A, bl, bu, xl, xu
diff --git a/test/runtests.jl b/test/runtests.jl
index 2b35495..35c9b95 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,6 +9,6 @@ const MatOI = MatrixOptInterface
 const MOI = MatOI.MOI
 # CUDA.allowscalar(false)
 
-include("PhaseOne/Artificial.jl")
-include("PhaseOne/Cplex.jl")
+# include("PhaseOne/Artificial.jl")
+# include("PhaseOne/Cplex.jl")
 include("netlib.jl")