diff --git a/config.ld b/config.ld
new file mode 100644
index 00000000..371884f1
--- /dev/null
+++ b/config.ld
@@ -0,0 +1,10 @@
+file = "modfiles"
+project = "Factory Planner"
+title = "Factory Planner Refarence"
+all = true
+format = "markdown"
+dir = "docs"
+
+custom_see_handler("^(https?://.*)$", function(url)
+    return url, url
+end)
\ No newline at end of file
diff --git a/modfiles/control.lua b/modfiles/control.lua
index 877af21f..b1d09a79 100644
--- a/modfiles/control.lua
+++ b/modfiles/control.lua
@@ -5,20 +5,10 @@ NTH_TICK_HANDLERS = {}
 GENERIC_HANDLERS = {}
 SEARCH_HANDLERS = {}
 
-require("util")  -- core.lualib
-table = require('__flib__.table')  -- replaces the lua table module
-
-require("data.init")
-require("data.data_util")
-
-require("ui.dialogs.main_dialog")
-require("ui.dialogs.modal_dialog")
-require("ui.ui_util")
-require("ui.event_handler")
-
 DEVMODE = true  -- Enables certain conveniences for development
 MARGIN_OF_ERROR = 1e-8  -- Margin of error for floating point calculations
 TIMESCALE_MAP = {[1] = "second", [60] = "minute", [3600] = "hour"}
+SOLVER_TYPE_MAP = {[1] = "traditional", [2] = "matrix", [3] = "interior_point"}
 SUBFACTORY_DELETION_DELAY = 15 * 60 * 60 -- ticks to deletion after subfactory trashing
 NEW = nil  -- global variable used to store new prototype data temporarily for migration
 
@@ -46,6 +36,17 @@ if DEVMODE then
     DEV_EXPORT_STRING = "eNq1lFFr2zAQx7+Lnu3gpE0Zfh0bFFYY22MpRpbPyRXJ8uRzwQR/951kmaypN1a3fYvv/rn76X86nYSxVfEErkPbiFxkm+2nzfVOJKLry1oqsg6hE/n9STTSACs4hcprT4KG1keQwHA05tHZJm21JBBjIggNdEpqztxkrLHkq/ki352tekW+ji0fQdHUpXWWrA/GcqA5xSVRpQqd6pE8G5pWY41QiZxcD8kzFG7r4FePDqpCGts3oUkFNTYcKQfWxXAy/8ivs8zT2rbQ8AR6Lqu07DzvDDsmLwl7JxvsTbq72n8I2nY1WRiFdfAeWCVontCZKtvskxAsLrsS+9G11lHq0wu9x7Xn0Xg4UmpR/+tAte6xWmX0br/a6Y5Amo+h2mb/eTUfkrhIxZy6nRZz/vxstV8mv+ixYK2tdR7hGwNcbuL8t5Dz4ArbIFq5ooqfhIN1/pDKyZqwOXiOduJnqiJ6NUXgD/IfU29WG6mOkfUSg6VgSs1l06hKd68BOUpXFRoN4+e11B0r7xhEPzfmNL6wOar+ZrQJ6ULFvTnL7uJZvAngFDQkDxAG7l1BfnlpKKI985YFsIVRxyEtE8R7sz3nv4bBj+fAz/m1H8L9jpZ+MS0N4rWv+IJFb7iJC9Xe87AP428/DGfh"
 end
 
+require("util")  -- core.lualib
+table = require('__flib__.table')  -- replaces the lua table module
+
+require("data.init")
+require("data.data_util")
+
+require("ui.dialogs.main_dialog")
+require("ui.dialogs.modal_dialog")
+require("ui.ui_util")
+require("ui.event_handler")
+
 
 -- ** UTIL **
 -- No better place for this too simple, yet too specific function anywhere else
diff --git a/modfiles/data/calculation/Matrix.lua b/modfiles/data/calculation/Matrix.lua
new file mode 100644
index 00000000..5128b84b
--- /dev/null
+++ b/modfiles/data/calculation/Matrix.lua
@@ -0,0 +1,366 @@
+local P, S = {}, {}
+local C = require("data.calculation.class")
+
+function P:__new(height, width)
+    self.height = height
+    self.width = width
+    for y = 1, height do
+        self[y] = {}
+    end
+end
+
+function S.new_vector(degree)
+    return S(degree, 1)
+end
+
+function S.list_to_vector(list, degree)
+    degree = degree or #list
+    local ret = S.new_vector(degree)
+    for y = 1, degree do
+        ret[y][1] = list[y] or 0
+    end
+    return ret
+end
+
+function S.diag(vector)
+    assert(vector.width == 1)
+    local size = vector.height
+    local ret = S(size, size):fill(0)
+    for i = 1, size do
+        ret[i][i] = vector[i][1]
+    end
+    return ret
+end
+
+function S.join(matrixes)
+    local heights, widths = {}, {}
+    for y, t in ipairs(matrixes) do
+        assert(#matrixes[1] == #t)
+        for x, v in ipairs(t) do
+            if S.is_matrix(v) then
+                if not heights[y] then
+                    heights[y] = v.height
+                else
+                    assert(heights[y] == v.height)
+                end
+                if not widths[x] then
+                    widths[x] = v.width
+                else
+                    assert(widths[x] == v.width)
+                end
+            else
+                assert(type(v) == "number")
+            end
+        end
+    end
+    
+    local total_height, total_width = 0, 0
+    for i = 1, #matrixes do
+        heights[i] = heights[i] or 1
+        total_height = total_height + heights[i]
+    end
+    for i = 1, #matrixes[1] do
+        widths[i] = widths[i] or 1
+        total_width = total_width + widths[i]
+    end
+    
+    local ret = S(total_height, total_width)
+    local y_offset = 0
+    for oy, t in ipairs(matrixes) do
+        local x_offset = 0
+        for ox, m in ipairs(t) do
+            if S.is_matrix(m) then
+                for y = 1, m.height do
+                    for x = 1, m.width do
+                        ret[y + y_offset][x + x_offset] = m[y][x]
+                    end
+                end
+            else
+                assert(m == 0 or widths[ox] == heights[oy])
+                for y = 1, heights[oy] do
+                    for x = 1, widths[ox] do
+                        if y == x then
+                            ret[y + y_offset][x + x_offset] = m
+                        else
+                            ret[y + y_offset][x + x_offset] = 0
+                        end
+                    end
+                end
+            end
+            x_offset = x_offset + widths[ox]
+        end
+        y_offset = y_offset + heights[oy]
+    end
+    return ret
+end
+
+function S.join_vector(vectors)
+    local matrixes = {}
+    for i, v in ipairs(vectors) do
+        matrixes[i] = {v}
+    end
+    return S.join(matrixes)
+end
+
+function S.is_matrix(value)
+    return C.class_type(value) == C.class_type(S)
+end
+
+function P:clone()
+    local height, width = self.height, self.width
+    local ret = S(height, width)
+    for y = 1, height do
+        for x = 1, width do
+            ret[y][x] = self[y][x]
+        end
+    end
+    return ret
+end
+
+function P:fill(value)
+    local height, width = self.height, self.width
+    for y = 1, height do
+        for x = 1, width do
+            self[y][x] = value
+        end
+    end
+    return self
+end
+
+function P.__add(op1, op2)
+    assert(S.is_matrix(op1) and S.is_matrix(op2))
+    assert(op1.height == op2.height and op1.width == op2.width)
+    local height, width = op1.height, op1.width
+    local ret = S(height, width)
+    for y = 1, height do
+        for x = 1, width do
+            ret[y][x] = op1[y][x] + op2[y][x]
+        end
+    end
+    return ret
+end
+
+function P.__sub(op1, op2)
+    assert(S.is_matrix(op1) and S.is_matrix(op2))
+    assert(op1.height == op2.height and op1.width == op2.width)
+    local height, width = op1.height, op1.width
+    local ret = S(height, width)
+    for y = 1, height do
+        for x = 1, width do
+            ret[y][x] = op1[y][x] - op2[y][x]
+        end
+    end
+    return ret
+end
+
+function P.__mul(op1, op2)
+    local function mul_scalar(m, s)
+        local height, width = m.height, m.width
+        local ret = S(height, width)
+        for y = 1, height do
+            for x = 1, width do
+                ret[y][x] = m[y][x] * s
+            end
+        end
+        return ret
+    end
+    
+    if type(op1) == "number" then
+        return mul_scalar(op2, op1)
+    elseif type(op2) == "number" then
+        return mul_scalar(op1, op2)
+    elseif S.is_matrix(op1) and S.is_matrix(op2) then
+        assert(op1.width == op2.height)
+        local l, height, width = op1.width, op1.height, op2.width
+        local ret = S(height, width)
+        for y = 1, height do
+            for x = 1, width do
+                local v = 0
+                for r = 1, l do
+                    v = v + op1[y][r] * op2[r][x]
+                end
+                ret[y][x] = v
+            end
+        end
+        return ret
+    else
+        assert()
+    end
+end
+
+function P.__div(op1, op2)
+    if type(op1) == "number" then
+        local height, width = op2.height, op2.width
+        local ret = S(height, width)
+        for y = 1, height do
+            for x = 1, width do
+                ret[y][x] = op1 / op2[y][x]
+            end
+        end
+        return ret
+    elseif type(op2) == "number" then
+        local height, width = op1.height, op1.width
+        local ret = S(height, width)
+        for y = 1, height do
+            for x = 1, width do
+                ret[y][x] = op1[y][x] / op2
+            end
+        end
+        return ret
+    else
+        assert()
+    end
+end
+
+function P:__unm()
+    local height, width = self.height, self.width
+    local ret = S(height, width)
+    for y = 1, height do
+        for x = 1, width do
+            ret[y][x] = -self[y][x]
+        end
+    end
+    return ret
+end
+
+function S.hadamard_product(op1, op2)
+    assert(S.is_matrix(op1) and S.is_matrix(op2))
+    assert(op1.height == op2.height and op1.width == op2.width)
+    local height, width = op1.height, op1.width
+    local ret = S(height, width)
+    for y = 1, height do
+        for x = 1, width do
+            ret[y][x] = op1[y][x] * op2[y][x]
+        end
+    end
+    return ret
+end
+
+function S.hadamard_power(matrix, scalar)
+    assert(S.is_matrix(matrix) and type(scalar) == "number")
+    local height, width = matrix.height, matrix.width
+    local ret = S(height, width)
+    for y = 1, height do
+        for x = 1, width do
+            local v = matrix[y][x]
+            ret[y][x] = v ^ scalar
+        end
+    end
+    return ret
+end
+
+function P:T()
+    local height, width = self.width, self.height
+    local ret = S(height, width)
+    for y = 1, height do
+        for x = 1, width do
+            ret[y][x] = self[x][y]
+        end
+    end
+    return ret
+end
+
+function P:sum()
+    local height, width = self.height, self.width
+    assert(width == 1)
+    local ret = 0
+    for i = 1, height do
+        ret = ret + self[i][1]
+    end
+    return ret
+end
+
+function P:euclidean_norm()
+    local height, width = self.height, self.width
+    assert(width == 1)
+    local ret = 0
+    for i = 1, height do
+        ret = ret + self[i][1] ^ 2
+    end
+    return math.sqrt(ret)
+end
+
+function P:submatrix(top, left, bottom, right)
+    local height, width = 1 + bottom - top, 1 + right - left
+    local ret = S(height, width)
+    for y = 1, height do
+        for x = 1, width do
+            ret[y][x] = self[top + y - 1][left + x - 1]
+        end
+    end
+    return ret
+end
+
+function P:insert_column(vector, x)
+    self.width = self.width + 1
+    x = x or self.width
+    if type(vector) == "number" then
+        local value = vector
+        for y = 1, self.height do
+            table.insert(self[y], x, value)
+        end
+    else
+        assert(S.is_matrix(vector) and vector.width == 1)
+        for y = 1, self.height do
+            table.insert(self[y], x, vector[y][1])
+        end
+    end
+    return self
+end
+
+function P:remove_column(x)
+    x = x or self.width
+    self.width = self.width - 1
+    local ret = {}
+    for y = 1, self.height do
+        ret[y] = table.remove(self[y], x)
+    end
+    return S.list_to_vector(ret)
+end
+
+function P:get(y, x)
+    return self[y][x]
+end
+
+function P:set(y, x, value)
+    self[y][x] = value
+end
+
+function P:iterate_row(y)
+    local x, width = 0, self.width
+    local function it()
+        x = x + 1
+        if x <= width then
+            return x, self[y][x]
+        else
+            return nil
+        end
+    end
+    return it
+end
+
+function P:row_swap(a, b)
+    self[a], self[b] = self[b], self[a]
+    return self
+end
+
+function P:row_mul(y, factor)
+    for x, v in ipairs(self[y]) do
+        self[y][x] = v * factor
+    end
+    return self
+end
+
+function P:row_trans(to, from, factor)
+    assert(to ~= from)
+    if factor == 0 then
+        return self
+    end
+    factor = factor or 1
+    for x, v in ipairs(self[from]) do
+        self[to][x] = self[to][x] + v * factor
+    end
+    return self
+end
+
+return C.class("Matrix", P, S)
diff --git a/modfiles/data/calculation/Problem.lua b/modfiles/data/calculation/Problem.lua
new file mode 100644
index 00000000..0a121708
--- /dev/null
+++ b/modfiles/data/calculation/Problem.lua
@@ -0,0 +1,247 @@
+--- Helper for generating linear programming problems.
+-- @classmod Problem
+-- @alias P
+-- @license MIT
+-- @author B_head
+
+local P, S = {}, {}
+local C = require("data.calculation.class")
+local Matrix = require("data.calculation.Matrix")
+local SparseMatrix = require("data.calculation.SparseMatrix")
+
+--- Constructor.
+-- @tparam string name Name of the problem.
+function P:__new(name)
+    self.name = name
+    self.primal = {}
+    self.primal_length = 0
+    self.dual = {}
+    self.dual_length = 0
+    self.subject_terms = {}
+end
+
+--- Add the variables to optimize and a term for the objective function.
+-- @tparam string key Term key.
+-- @tparam number cost Coefficient of the term.
+-- @tparam bool is_result If ture, variables are included in the output of the solution.
+function P:add_objective_term(key, cost, is_result)
+    assert(not self.primal[key])
+    self.primal_length = self.primal_length + 1
+    self.primal[key] = {
+        key = key,
+        index = self.primal_length,
+        value = cost or 1,
+        is_result = is_result or false,
+    }
+end
+
+--- Is there an objective term that corresponds to the key?
+-- @tparam string key Term key.
+-- @treturn bool True if exists.
+function P:is_exist_objective(key)
+    return self.primal[key] ~= nil
+end
+
+--- Add to the coefficients of the terms of the objective function.
+-- @tparam string key Term key.
+-- @tparam number cost Value to be added.
+function P:add_objective_cost(key, cost)
+    assert(self.primal[key])
+    self.primal[key].value = self.primal[key].value + cost
+end
+
+--- Add equality constraints.
+-- @tparam string key Constraint key.
+-- @tparam number target The value that needs to match the result of the constraint expression.
+function P:add_eq_constraint(key, target)
+    assert(not self.dual[key])
+    self.dual_length = self.dual_length + 1
+    self.dual[key] = {
+        key = key,
+        index = self.dual_length,
+        value = target or 0,
+    }
+end
+
+--- Add inequality constraints of equal or less.
+-- @tparam string key Constraint key.
+-- @tparam number limit The upper bound on the result of the constraint equation.
+function P:add_le_constraint(key, limit)
+    local slack_key = "<slack>" .. key
+    self:add_eq_constraint(key, limit)
+    self:add_objective_term(slack_key, 0, false)
+    self:add_constraint_term(slack_key, {
+        [key] = 1,
+    })
+end
+
+--- Add inequality constraints of equal or greater.
+-- @tparam string key Constraint key.
+-- @tparam number limit The lower bound on the result of the constraint equation.
+function P:add_ge_constraint(key, limit)
+    local slack_key = "<slack>" .. key
+    self:add_eq_constraint(key, limit)
+    self:add_objective_term(slack_key, 0, false)
+    self:add_constraint_term(slack_key, {
+        [key] = -1,
+    })
+end
+
+--- Is there an constraint that corresponds to the key?
+-- @tparam string key Constraint key.
+-- @treturn bool True if exists.
+function P:is_exist_constraint(key)
+    return self.dual[key] ~= nil
+end
+
+--- Add a term for the constraint equation.
+-- @tparam string objective_key Objective term key.
+-- @tparam {[string]=number,...} subject_map Mapping of constraint keys to term coefficients.
+function P:add_constraint_term(objective_key, subject_map)
+    local term = self.subject_terms[objective_key] or {}
+    for k, v in pairs(subject_map) do
+        assert(not term[k])
+        term[k] = v
+    end
+    self.subject_terms[objective_key] = term
+end
+
+--- Make a vector of the coefficient of the primal problem.
+-- @treturn Matrix Coefficients in vector form.
+function P:make_primal_coefficients()
+    local ret = Matrix.new_vector(self.primal_length)
+    for _, v in pairs(self.primal) do
+        ret[v.index][1] = v.value
+    end
+    return ret
+end
+
+--- Make a vector of the coefficient of the dual problem.
+-- @treturn Matrix Coefficients in vector form.
+function P:make_dual_coefficients()
+    local ret = Matrix.new_vector(self.dual_length)
+    for _, v in pairs(self.dual) do
+        ret[v.index][1] = v.value
+    end
+    return ret
+end
+
+--- Make a sparse matrix of constraint equations.
+-- @treturn SparseMatrix Constraint equations in matrix form.
+function P:make_subject_sparse_matrix()
+    local ret = SparseMatrix(self.dual_length, self.primal_length)
+    for p, t in pairs(self.subject_terms) do
+        if self.primal[p] then
+            local x = self.primal[p].index
+            for d, v in pairs(t) do
+                if self.dual[d] then
+                    local y = self.dual[d].index
+                    ret:set(y, x, v)
+                end
+            end
+        end
+    end
+    return ret
+end
+
+--- Make a matrix of constraint equations.
+-- @treturn Matrix Constraint equations in matrix form.
+function P:make_subject_matrix()
+    return self:make_subject_sparse_matrix():to_matrix()
+end
+
+--- Make find variables for primal problem.
+-- @tparam table prev_raw_solution The value returned by @{pack_pdip_variables}.
+-- @treturn Matrix Variables in vector form.
+function P:make_primal_find_variables(prev_raw_solution)
+    local prev_x = prev_raw_solution and prev_raw_solution.x or {}
+    local ret = Matrix.new_vector(self.primal_length)
+    for k, v in pairs(self.primal) do
+        ret[v.index][1] = prev_x[k] or 1
+    end
+    return ret
+end
+
+--- Make find variables for dual problem.
+-- @tparam table prev_raw_solution The value returned by @{pack_pdip_variables}.
+-- @treturn Matrix Variables in vector form.
+function P:make_dual_find_variables(prev_raw_solution)
+    local prev_y = prev_raw_solution and prev_raw_solution.y or {}
+    local ret = Matrix.new_vector(self.dual_length)
+    for k, v in pairs(self.dual) do
+        ret[v.index][1] = prev_y[k] or 0
+    end
+    return ret
+end
+
+--- Make slack variables for dual problem.
+-- @tparam table prev_raw_solution The value returned by @{pack_pdip_variables}.
+-- @treturn Matrix Variables in vector form.
+function P:make_dual_slack_variables(prev_raw_solution)
+    local prev_s = prev_raw_solution and prev_raw_solution.s or {}
+    local ret = Matrix.new_vector(self.primal_length)
+    for k, v in pairs(self.primal) do
+        ret[v.index][1] = prev_s[k] or math.max(1, v.value)
+    end
+    return ret
+end
+
+--- Store the value of variables in a plain table.
+-- @tparam Matrix x Find variables for primal problem.
+-- @tparam Matrix y Find variables for dual problem.
+-- @tparam Matrix s Slack variables for dual problem.
+-- @treturn table Packed table.
+function P:pack_pdip_variables(x, y, s)
+    local ret_x, ret_s = {}, {}
+    for k, v in pairs(self.primal) do
+        ret_x[k] = x[v.index][1]
+        ret_s[k] = s[v.index][1]
+    end
+
+    local ret_y = {}
+    for k, v in pairs(self.dual) do
+        ret_y[k] = y[v.index][1]
+    end
+    
+    return {
+        x = ret_x,
+        y = ret_y,
+        s = ret_s,
+    }
+end
+
+--- Remove unnecessary variables from the solution of the problem.
+-- @tparam Matrix vector Solution to the primal problem.
+-- @treturn table Filtered solution.
+function P:filter_solution_to_result(vector)
+    local ret = {}
+    for k, v in pairs(self.primal) do
+        if v.is_result then
+            ret[k] = vector[v.index][1]
+        end
+    end
+    return ret
+end
+
+--- Put primal variables in readable format.
+-- @tparam Matrix vector Variables to the primal problem.
+-- @treturn string Readable text.
+function P:dump_primal(vector)
+    local ret = ""
+    for k, v in pairs(self.primal) do
+        ret = ret .. string.format("%q = %f\n", k, vector[v.index][1])
+    end
+    return ret
+end
+
+--- Put dual variables in readable format.
+-- @treturn string Readable text.
+function P:dump_dual(vector)
+    local ret = ""
+    for k, v in pairs(self.dual) do
+        ret = ret .. string.format("%q = %f\n", k, vector[v.index][1])
+    end
+    return ret
+end
+
+return C.class("Problem", P, S)
\ No newline at end of file
diff --git a/modfiles/data/calculation/SparseMatrix.lua b/modfiles/data/calculation/SparseMatrix.lua
new file mode 100644
index 00000000..c133a463
--- /dev/null
+++ b/modfiles/data/calculation/SparseMatrix.lua
@@ -0,0 +1,360 @@
+local P, S = {}, {}
+local C = require("data.calculation.class")
+local Matrix = require("data.calculation.Matrix")
+
+function P:__new(height, width)
+    self.height = height
+    self.width = width
+    self.values = {}
+    self.indexes = {}
+    for y = 1, height do
+        self.values[y] = {}
+        self.indexes[y] = {}
+    end
+end
+
+function S.diag(vector)
+    assert(vector.width == 1)
+    local size = vector.height
+    local ret = S(size, size)
+    for i = 1, size do
+        ret:set(i, i, vector[i][1])
+    end
+    return ret
+end
+
+function S.join(matrixes)
+    local heights, widths = {}, {}
+    for y, t in ipairs(matrixes) do
+        assert(#matrixes[1] == #t)
+        for x, v in ipairs(t) do
+            if S.is_sparse_matrix(v) then
+                if not heights[y] then
+                    heights[y] = v.height
+                else
+                    assert(heights[y] == v.height)
+                end
+                if not widths[x] then
+                    widths[x] = v.width
+                else
+                    assert(widths[x] == v.width)
+                end
+            else
+                assert(type(v) == "number")
+            end
+        end
+    end
+    
+    local total_height, total_width = 0, 0
+    for i = 1, #matrixes do
+        heights[i] = heights[i] or 1
+        total_height = total_height + heights[i]
+    end
+    for i = 1, #matrixes[1] do
+        widths[i] = widths[i] or 1
+        total_width = total_width + widths[i]
+    end
+    
+    local ret = S(total_height, total_width)
+    local y_offset = 0
+    for oy, t in ipairs(matrixes) do
+        local x_offset = 0
+        for ox, m in ipairs(t) do
+            if S.is_sparse_matrix(m) then
+                for y = 1, m.height do
+                    local values = m.values[y]
+                    for rx, x in ipairs(m.indexes[y]) do
+                        table.insert(ret.values[y + y_offset], values[rx])
+                        table.insert(ret.indexes[y + y_offset], x + x_offset)
+                    end
+                end
+            elseif m ~= 0 then
+                assert(widths[ox] == heights[oy])
+                local size = heights[oy]
+                for i = 1, size do
+                    table.insert(ret.values[i + y_offset], m)
+                    table.insert(ret.indexes[i + y_offset], i + x_offset)
+                end
+            end
+            x_offset = x_offset + widths[ox]
+        end
+        y_offset = y_offset + heights[oy]
+    end
+    return ret
+end
+
+function S.is_sparse_matrix(value)
+    return C.class_type(value) == C.class_type(S)
+end
+
+function P:to_matrix()
+    local height, width = self.height, self.width
+    local ret = Matrix(height, width):fill(0)
+    for y = 1, height do
+        local values = self.values[y]
+        for rx, v in ipairs(self.indexes[y]) do
+            ret[y][v] = values[rx]
+        end
+    end
+    return ret
+end
+
+function P:clone()
+    local height, width = self.height, self.width
+    local ret = S(height, width)
+    for y = 1, height do
+        for rx, v in ipairs(self.values[y]) do
+            ret.values[y][rx] = v
+        end
+        for rx, v in ipairs(self.indexes[y]) do
+            ret.indexes[y][rx] = v
+        end
+    end
+    return ret
+end
+
+function P.__mul(op1, op2)
+    local function mul_scalar(m, s)
+        local height, width = m.height, m.width
+        local ret = S(height, width)
+        for y = 1, height do
+            for rx, v in ipairs(m.values[y]) do
+                ret.values[y][rx] = v * s
+            end
+            for rx, v in ipairs(m.indexes[y]) do
+                ret.indexes[y][rx] = v
+            end
+        end
+        return ret
+    end
+    
+    if type(op1) == "number" then
+        return mul_scalar(op2, op1)
+    elseif type(op2) == "number" then
+        return mul_scalar(op1, op2)
+    elseif S.is_sparse_matrix(op1) and Matrix.is_matrix(op2) then
+        assert(op1.width == op2.height)
+        local height, width = op1.height, op2.width
+        local ret = Matrix(height, width)
+        for y = 1, height do
+            for x = 1, width do
+                local op1_row_indexes, op1_row_values = op1.indexes[y], op1.values[y]
+                local v = 0
+                for rx, r in ipairs(op1_row_indexes) do
+                    v = v + op1_row_values[rx] * op2[r][x]
+                end
+                ret[y][x] = v
+            end
+        end
+        return ret
+    elseif Matrix.is_matrix(op1) and S.is_sparse_matrix(op2) then
+        assert(op1.width == op2.height)
+        local height, width = op1.height, op2.width
+        local op2_t = op2:T()
+        local ret = Matrix(height, width)
+        for y = 1, height do
+            for x = 1, width do
+                local op2_column_indexes, op2_column_values = op2_t.indexes[x], op2_t.values[x]
+                local v = 0
+                for ry, r in ipairs(op2_column_indexes) do
+                    v = v + op1[y][r] * op2_column_values[ry]
+                end
+                ret[y][x] = v
+            end
+        end
+        return ret
+    elseif S.is_sparse_matrix(op1) and S.is_sparse_matrix(op2) then
+        assert(op1.width == op2.height)
+        local height, width = op1.height, op2.width
+        local op2_t = op2:T()
+        local ret = S(height, width)
+        for y = 1, height do
+            local ret_indexes, ret_values = ret.indexes[y], ret.values[y]
+            for x = 1, width do
+                local op1_row_indexes, op1_row_values = op1.indexes[y], op1.values[y]
+                local op2_column_indexes, op2_column_values = op2_t.indexes[x], op2_t.values[x]
+                local ry, rx = 1, 1
+                local function get_x()
+                    return op1_row_indexes[rx] or math.huge, op2_column_indexes[ry] or math.huge
+                end
+
+                local v = 0
+                local op1_r, op2_r = get_x()
+                while not (op1_r == math.huge and op2_r == math.huge) do
+                    if op1_r < op2_r then
+                        rx = rx + 1
+                    elseif op1_r > op2_r then
+                        ry = ry + 1
+                    else -- op1_r == op2_r
+                        v = v + op1_row_values[rx] * op2_column_values[ry]
+                        rx = rx + 1
+                        ry = ry + 1
+                    end
+                    op1_r, op2_r = get_x()
+                end
+                if v ~= 0 then
+                    table.insert(ret_values, v)
+                    table.insert(ret_indexes, x)
+                end
+            end
+        end
+        return ret
+    else
+        assert()
+    end
+end
+
+function P:__unm()
+    return self:__mul(-1)
+end
+
+function P:T()
+    local height, width = self.width, self.height
+    local ret = S(height, width)
+    for x = 1, width do
+        local values = self.values[x]
+        for rx, y in ipairs(self.indexes[x]) do
+            table.insert(ret.values[y], values[rx])
+            table.insert(ret.indexes[y], x)
+        end
+    end
+    return ret
+end
+
+function P:get_raw_index(y, x)
+    local indexes = self.indexes[y]
+    if x > self.width then
+        return #indexes + 1
+    end
+    for rx, v in ipairs(indexes) do
+        if x <= v then
+            return rx, x == v
+        end
+    end
+    return #indexes + 1
+end
+
+function P:insert_column(vector, x)
+    x = x or self.width + 1
+    if type(vector) == "number" then
+        local value = vector
+        for y = 1, self.height do
+            self:set(y, x, value)
+        end
+    else
+        assert(Matrix.is_matrix(vector) and vector.width == 1)
+        for y = 1, self.height do
+            self:set(y, x, vector[y][1])
+        end
+    end
+    self.width = self.width + 1
+    return self
+end
+
+function P:remove_column(x)
+    x = x or self.width
+    local ret = {}
+    for y = 1, self.height do
+        local rx, e = self:get_raw_index(y, x)
+        if e then
+            local values, indexes = self.values[y], self.indexes[y]
+            ret[y] = values[rx]
+            table.remove(values, rx)
+            table.remove(indexes, rx)
+        else
+            ret[y] = 0
+        end
+    end
+    self.width = self.width - 1
+    return Matrix.list_to_vector(ret)
+end
+
+function P:get(y, x)
+    local rx, e = self:get_raw_index(y, x)
+    if e then
+        return self.values[y][rx]
+    else
+        return 0
+    end
+end
+
+function P:set(y, x, value)
+    local rx, e = self:get_raw_index(y, x)
+    if e then
+        if value ~= 0 then
+            self.values[y][rx] = value
+        else
+            table.remove(self.values[y], rx)
+            table.remove(self.indexes[y], rx)
+        end
+    elseif value ~= 0 then
+        table.insert(self.values[y], rx, value)
+        table.insert(self.indexes[y], rx, x)
+    end
+end
+
+function P:iterate_row(y)
+    local values, indexes = self.values[y], self.indexes[y]
+    local rx = 0
+    local function it()
+        rx = rx + 1
+        return indexes[rx], values[rx]
+    end
+    return it
+end
+
+function P:row_swap(a, b)
+    self.values[a], self.values[b] = self.values[b], self.values[a]
+    self.indexes[a], self.indexes[b] = self.indexes[b], self.indexes[a]
+    return self
+end
+
+function P:row_mul(y, factor)
+    local values = self.values[y]
+    for rx, v in ipairs(values) do
+        values[rx] = v * factor
+    end
+    return self
+end
+
+function P:row_trans(to, from, factor)
+    assert(to ~= from)
+    if factor == 0 then
+        return self
+    end
+    factor = factor or 1
+    local to_values, to_indexes, to_rx = self.values[to], self.indexes[to], 1
+    local from_values, from_indexes, from_rx = self.values[from], self.indexes[from], 1
+    local new_values, new_indexes = {}, {}
+    local function get_x()
+        return to_indexes[to_rx] or math.huge, from_indexes[from_rx] or math.huge
+    end
+
+    local to_x, from_x = get_x()
+    while not (to_x == math.huge and from_x == math.huge) do
+        if to_x < from_x then
+            local v = to_values[to_rx]
+            table.insert(new_values, v)
+            table.insert(new_indexes, to_x)
+            to_rx = to_rx + 1
+        elseif to_x > from_x then
+            local v = from_values[from_rx] * factor
+            table.insert(new_values, v)
+            table.insert(new_indexes, from_x)
+            from_rx = from_rx + 1
+        else -- to_x == from_x
+            local v = to_values[to_rx] + from_values[from_rx] * factor
+            if v ~= 0 then
+                table.insert(new_values, v)
+                table.insert(new_indexes, to_x)
+            end
+            to_rx = to_rx + 1
+            from_rx = from_rx + 1
+        end
+        to_x, from_x = get_x()
+    end
+    self.values[to], self.indexes[to] = new_values, new_indexes
+    return self
+end
+
+return C.class("SparseMatrix", P, S)
\ No newline at end of file
diff --git a/modfiles/data/calculation/class.lua b/modfiles/data/calculation/class.lua
new file mode 100644
index 00000000..a891a5ba
--- /dev/null
+++ b/modfiles/data/calculation/class.lua
@@ -0,0 +1,95 @@
+--- Helper for generate OOP-style data structures.
+-- @module class
+-- @license MIT
+-- @author B_head
+
+local M = {}
+
+local function noop()
+    -- no operation.
+end
+
+local function create_metatable(name, prototype, extend_class)
+    local super_metatable = getmetatable(extend_class) or {}
+    local super_prototype = super_metatable.__prototype
+    setmetatable(prototype, {
+        __index = super_prototype -- Constructing prototype chains.
+    })
+    local ret = {
+        __new = noop
+    }
+    for k, v in pairs(super_metatable) do
+        ret[k] = v
+    end
+    ret.__type = name
+    ret.__prototype = prototype
+    ret.__super_prototype = super_prototype
+    ret.__extend = extend_class
+    ret.__index = prototype -- Assign a prototype to instances.
+    for k, v in pairs(prototype) do
+        if k:sub(1, 2) == "__" then
+            ret[k] = v
+        end
+    end
+    return ret
+end
+
+local function create_instance(class_object, ...)
+    local mt = getmetatable(class_object)
+    local ret = {}
+    setmetatable(ret, mt)
+    mt.__new(ret, ...)
+    return ret
+end
+
+--- Create class object.
+-- @tparam string name Name of the class type.
+-- @tparam table prototype A table that defines methods, meta-methods, and constants.
+-- @tparam table static A table that defines static functions.
+-- @param extend_class Class object to inherit from.
+-- @return Class object.
+function M.class(name, prototype, static, extend_class)
+    static = static or {}
+    setmetatable(static, {
+        __call = create_instance,
+        -- Overrides the metatable returned by getmetatable(class_object).
+        __metatable = create_metatable(name, prototype, extend_class), 
+    })
+    return static -- Return as class_object.
+end
+
+--- Return name of the class type.
+-- @param value Class object.
+-- @return Class name.
+function M.class_type(value)
+    local mt = getmetatable(value)
+    return mt and mt.__type
+end
+
+--- Return a prototype table.
+-- @param value Class object.
+-- @return Prototype table.
+function M.prototype(value)
+    local mt = getmetatable(value)
+    return mt and mt.__prototype
+end
+
+--- Return a prototype table of the superclass.
+-- @param value Class object.
+-- @return Prototype table.
+function M.super(value)
+    local mt = getmetatable(value)
+    return mt and mt.__super_prototype
+end
+
+--- Restore methods, meta-methods, and constants in the instance table.
+-- @tparam table plain_table An instance table to restore.
+-- @param class_object Class object that defines methods, meta-methods, and constants.
+-- @return Instance table.
+function M.resetup(plain_table, class_object)
+    local mt = getmetatable(class_object)
+    setmetatable(plain_table, mt)
+    return plain_table
+end
+
+return M
\ No newline at end of file
diff --git a/modfiles/data/calculation/interface.lua b/modfiles/data/calculation/interface.lua
index 4ee096c6..5f02bf4e 100644
--- a/modfiles/data/calculation/interface.lua
+++ b/modfiles/data/calculation/interface.lua
@@ -177,6 +177,8 @@ local function update_ingredient_satisfaction(floor, product_class)
     end
 end
 
+local solver_util = require("data.calculation.solver_util")
+local linear_optimization_solver = require("data.calculation.linear_optimization_solver")
 
 -- ** TOP LEVEL **
 -- Updates the whole subfactory calculations from top to bottom
@@ -187,10 +189,10 @@ function calculation.update(player, subfactory)
         player_table.active_subfactory = subfactory
 
         local subfactory_data = calculation.interface.generate_subfactory_data(player, subfactory)
-
-        if subfactory.matrix_free_items ~= nil then  -- meaning the matrix solver is active
+        if subfactory.solver_type == "matrix" then
             local matrix_metadata = matrix_solver.get_matrix_solver_metadata(subfactory_data)
 
+            subfactory.matrix_free_items = subfactory.matrix_free_items or {}
             if matrix_metadata.num_cols > matrix_metadata.num_rows and #subfactory.matrix_free_items > 0 then
                 subfactory.matrix_free_items = {}
                 subfactory_data = calculation.interface.generate_subfactory_data(player, subfactory)
@@ -214,8 +216,18 @@ function calculation.update(player, subfactory)
             else  -- reset top level items
                 set_blank_subfactory(player, subfactory)
             end
-        else
+        elseif subfactory.solver_type == "interior_point" then
+            local normalized_top_floor = solver_util.normalize(subfactory_data.top_floor)
+            local flat_recipe_lines = solver_util.to_flat_recipe_lines(normalized_top_floor)
+            local normalized_references = solver_util.normalize_references(subfactory_data.top_level_products, subfactory_data.timescale)
+            local problem = linear_optimization_solver.create_problem(subfactory_data.name, flat_recipe_lines, normalized_references)
+            local machine_counts, raw_solution = linear_optimization_solver.primal_dual_interior_point(problem, subfactory.prev_raw_solution)
+            solver_util.feedback(machine_counts, subfactory_data.player_index, subfactory_data.timescale, normalized_top_floor)
+            subfactory.prev_raw_solution = raw_solution
+        elseif subfactory.solver_type == "traditional" then
             sequential_solver.update_subfactory(subfactory_data)
+        else
+            assert(false, "Undefined solver_type = " .. subfactory.solver_type)
         end
 
         player_table.active_subfactory = nil  -- reset after calculations have been carried out
@@ -233,6 +245,8 @@ end
 function calculation.interface.generate_subfactory_data(player, subfactory)
     local subfactory_data = {
         player_index = player.index,
+        name = subfactory.name,
+        timescale = subfactory.timescale,
         top_level_products = {},
         top_floor = nil,
         matrix_free_items = subfactory.matrix_free_items
diff --git a/modfiles/data/calculation/linear_optimization_solver.lua b/modfiles/data/calculation/linear_optimization_solver.lua
new file mode 100644
index 00000000..33a4e6d3
--- /dev/null
+++ b/modfiles/data/calculation/linear_optimization_solver.lua
@@ -0,0 +1,539 @@
+--- Create and solve linear programming problems (a.k.a linear optimization).
+-- @module linear_optimization_solver
+-- @license MIT
+-- @author B_head
+
+local M = {}
+local Problem = require("data.calculation.Problem")
+local Matrix = require("data.calculation.Matrix")
+local SparseMatrix = require("data.calculation.SparseMatrix")
+
+local machine_count_penalty = 2 ^ 0
+local shortage_penalty = 2 ^ 25
+local surplusage_penalty = 2 ^ 15
+local products_priority_penalty = 2 ^ 5
+local ingredients_priority_penalty = 2 ^ 10
+
+local function get_include_items(flat_recipe_lines, normalized_references)
+    local set = {}
+    local function add_set(key, type)
+        if not set[key] then
+            set[key] = {
+                id = key,
+                product = false,
+                ingredient = false,
+                reference = false,
+            }
+        end
+        set[key][type] = true
+    end
+
+    for _, l in pairs(flat_recipe_lines) do
+        for k, _ in pairs(l.products) do
+            add_set(k, "product")
+        end
+        for k, _ in pairs(l.ingredients) do
+            add_set(k, "ingredient")
+        end
+    end
+    for k, _ in pairs(normalized_references) do
+        add_set(k, "reference")
+    end
+
+    return set
+end
+
+local function create_item_flow_graph(flat_recipe_lines)
+    local ret = {}
+    local function add(a, type, b, ratio)
+        if not ret[a] then
+            ret[a] = {
+                from = {},
+                to = {},
+                visited = false,
+                cycled = false,
+            }
+        end
+        table.insert(ret[a][type], {id=b, ratio=ratio})
+    end
+
+    for _, l in pairs(flat_recipe_lines) do
+        for _, a in pairs(l.products) do
+            for _, b in pairs(l.ingredients) do
+                local ratio = b.amount_per_machine_by_second / a.amount_per_machine_by_second
+                add(a.normalized_id, "to", b.normalized_id, ratio)
+            end
+        end
+        for _, a in pairs(l.ingredients) do
+            for _, b in pairs(l.products) do
+                local ratio = b.amount_per_machine_by_second / a.amount_per_machine_by_second
+                add(a.normalized_id, "from", b.normalized_id, ratio)
+            end
+        end
+    end
+    return ret
+end
+
+local function detect_cycle_dilemma_impl(item_flow_graph, id, path)
+    local current = item_flow_graph[id]
+    if current.visited then
+        local included = false
+        for _, path_id in ipairs(path) do
+            if path_id == id then
+                included = true
+            end
+            if included then
+                item_flow_graph[path_id].cycled = true
+            end
+        end
+        return
+    end
+
+    current.visited = true
+    table.insert(path, id)
+    for _, n in ipairs(current.to) do
+        detect_cycle_dilemma_impl(item_flow_graph, n.id, path)
+    end
+    table.remove(path)
+end
+
+local function detect_cycle_dilemma(flat_recipe_lines)
+    local item_flow_graph = create_item_flow_graph(flat_recipe_lines)
+    local path = {}
+    for id, _ in pairs(item_flow_graph) do
+        if not item_flow_graph[id].visited then
+            detect_cycle_dilemma_impl(item_flow_graph, id, path)
+        end
+    end
+
+    local ret = {}
+    for id, v in pairs(item_flow_graph) do
+        ret[id] = {product=v.cycled, ingredient=v.cycled}
+    end
+    return ret
+end
+
+--- Create linear programming problems.
+-- @tparam string problem_name Problem name.
+-- @param flat_recipe_lines List returned by @{solver_util.to_flat_recipe_lines}.
+-- @param normalized_references List returned by @{solver_util.normalize_references}.
+-- @treturn Problem Created problem object.
+function M.create_problem(problem_name, flat_recipe_lines, normalized_references)
+    local function add_item_factor(subject_map, name, factor)
+        subject_map["balance|" .. name] = factor
+        if factor > 0 then
+            subject_map["product_reference|" .. name] = factor
+        elseif factor < 0 then
+            subject_map["ingredient_reference|" .. name] = -factor
+        end
+    end
+
+    local problem = Problem(problem_name)
+    local need_slack = detect_cycle_dilemma(flat_recipe_lines)
+    for recipe_id, v in pairs(flat_recipe_lines) do
+        problem:add_objective_term(recipe_id, machine_count_penalty, true)
+        local subject_map = {}
+        local is_maximum_limit = false
+
+        if v.maximum_machine_count then
+            local key = "maximum|" .. recipe_id
+            problem:add_le_constraint(key, v.maximum_machine_count)
+            subject_map[key] = 1
+            is_maximum_limit = true
+        end
+        if v.minimum_machine_count then
+            local key = "minimum|" .. recipe_id
+            problem:add_ge_constraint(key, v.minimum_machine_count)
+            subject_map[key] = 1
+        end
+
+        for item_id, u in pairs(v.products) do
+            local amount = u.amount_per_machine_by_second
+            add_item_factor(subject_map, item_id, amount)
+            if is_maximum_limit then
+                need_slack[item_id].product = true
+            end
+
+            if #u.neighbor_recipe_lines >= 2 then
+                local balance_key = string.format("products_priority_balance|%s:%s", recipe_id, item_id)
+                problem:add_eq_constraint(balance_key, 0)
+                subject_map[balance_key] = amount
+                for priority, neighbor in ipairs(u.neighbor_recipe_lines) do
+                    local transfer_key = string.format("transfer|%s=>%s:%s", recipe_id, neighbor.normalized_id, item_id)
+                    local penalty = (priority - 1) * products_priority_penalty
+                    if problem:is_exist_objective(transfer_key) then
+                        problem:add_objective_cost(transfer_key, penalty)
+                    else
+                        problem:add_objective_term(transfer_key, penalty)
+                    end
+                    problem:add_constraint_term(transfer_key, {
+                        [balance_key] = -1
+                    })
+                end
+                local implicit_key = string.format("implicit_transfer|%s=>(?):%s", recipe_id, item_id)
+                local penalty = #u.neighbor_recipe_lines * products_priority_penalty
+                problem:add_objective_term(implicit_key, penalty)
+                problem:add_constraint_term(implicit_key, {
+                    [balance_key] = -1
+                })
+            end 
+        end
+
+        for item_id, u in pairs(v.ingredients) do
+            local amount = u.amount_per_machine_by_second
+            add_item_factor(subject_map, item_id, -amount)
+            if is_maximum_limit then
+                need_slack[item_id].ingredient = true
+            end
+
+            if #u.neighbor_recipe_lines >= 2 then
+                local balance_key = string.format("ingredients_priority_balance|%s:%s", recipe_id, item_id)
+                problem:add_eq_constraint(balance_key, 0)
+                subject_map[balance_key] = -amount
+                for priority, neighbor in ipairs(u.neighbor_recipe_lines) do
+                    local transfer_key = string.format("transfer|%s=>%s:%s", neighbor.normalized_id, recipe_id, item_id)
+                    local penalty = (priority - 1) * ingredients_priority_penalty
+                    if problem:is_exist_objective(transfer_key) then
+                        problem:add_objective_cost(transfer_key, penalty)
+                    else
+                        problem:add_objective_term(transfer_key, penalty)
+                    end
+                    problem:add_constraint_term(transfer_key, {
+                        [balance_key] = 1
+                    })
+                end
+                local implicit_key = string.format("implicit_transfer|(?)=>%s:%s", recipe_id, item_id)
+                local penalty = #u.neighbor_recipe_lines * ingredients_priority_penalty
+                problem:add_objective_term(implicit_key, penalty)
+                problem:add_constraint_term(implicit_key, {
+                    [balance_key] = 1
+                })
+            end 
+        end
+
+        problem:add_constraint_term(recipe_id, subject_map)
+    end
+    
+    local items = get_include_items(flat_recipe_lines, normalized_references)
+    for item_id, v in pairs(items) do
+        if v.product and v.ingredient then
+            problem:add_eq_constraint("balance|" .. item_id, 0)
+        end
+        if v.ingredient and need_slack[item_id].ingredient then
+            local key = "implicit_ingredient|" .. item_id
+            local penalty = surplusage_penalty
+            problem:add_objective_term(key, penalty)
+            local subject_map = {}
+            add_item_factor(subject_map, item_id, -1)
+            problem:add_constraint_term(key, subject_map)
+        end
+        if v.product and need_slack[item_id].product then
+            local key = "implicit_product|" .. item_id
+            local penalty = shortage_penalty
+            problem:add_objective_term(key, penalty)
+            local subject_map = {}
+            add_item_factor(subject_map, item_id, 1)
+            problem:add_constraint_term(key, subject_map)
+        end
+        if v.reference then
+            local r = normalized_references[item_id]
+            if v.product then
+                problem:add_ge_constraint("product_reference|" .. item_id, r.amount_per_second)
+            end
+            if v.ingredient then
+                problem:add_ge_constraint("ingredient_reference|" .. item_id, r.amount_per_second)
+            end
+        end
+    end
+
+    return problem
+end
+
+local debug_print = log
+local had, had_pow, diag = Matrix.hadamard_product, Matrix.hadamard_power, SparseMatrix.diag
+local tolerance = MARGIN_OF_ERROR
+local step_limit = 1 - (2 ^ -20)
+local machine_epsilon = (2 ^ -52)
+local tiny_number = math.sqrt(2 ^ -970)
+local iterate_limit = 200
+
+local function force_variables_constraint(variables)
+    local height = variables.height
+    for y = 1, height do
+        variables[y][1] = math.max(tiny_number, variables[y][1])
+    end
+end
+
+local function sigmoid(value, min, max)
+    min = min or 0
+    max = max or 1
+    return (max - min) / (1 + math.exp(-value)) + min
+end
+
+local function get_max_step(v, dir)
+    local height = v.height
+    local ret = 1
+    for y = 1, height do
+        local a, b = v[y][1], dir[y][1]
+        if b < 0 then
+            ret = math.min(ret, step_limit * (a / -b))
+        end
+    end
+    return ret
+end
+
+--- Solve linear programming problems.
+-- @see http://www.cas.mcmaster.ca/~cs777/presentations/NumericalIssue.pdf
+-- @tparam Problem problem Problems to solve.
+-- @tparam table prev_raw_solution The value returned by @{Problem:pack_pdip_variables}.
+-- @treturn {[string]=number,...} Solution to problem.
+-- @treturn table Packed table of raw solution.
+function M.primal_dual_interior_point(problem, prev_raw_solution)
+    local A = problem:make_subject_sparse_matrix()
+    local AT = A:T()
+    local b = problem:make_dual_coefficients()
+    local c = problem:make_primal_coefficients()
+    local p_degree = problem.primal_length
+    local d_degree = problem.dual_length
+    local x = problem:make_primal_find_variables(prev_raw_solution)
+    local y = problem:make_dual_find_variables(prev_raw_solution)
+    local s = problem:make_dual_slack_variables(prev_raw_solution)
+
+    debug_print(string.format("-- solve %s --", problem.name))
+    for i = 0, iterate_limit do
+        force_variables_constraint(x)
+        force_variables_constraint(s)
+
+        local dual = AT * y + s - c
+        local primal = A * x - b
+        local duality_gap = had(x, s)
+
+        local d_sat = (d_degree == 0) and 0 or dual:euclidean_norm() / d_degree
+        local p_sat = (p_degree == 0) and 0 or primal:euclidean_norm() / p_degree
+        local dg_sat = (p_degree == 0) and 0 or duality_gap:euclidean_norm() / p_degree
+
+        debug_print(string.format(
+            "i = %i, primal = %f, dual = %f, duality_gap = %f", 
+            i, p_sat, d_sat, dg_sat
+        ))
+        if math.max(d_sat, p_sat, dg_sat) <= tolerance then
+            break
+        end
+
+        local fvg = M.create_default_flee_value_generator(y)
+        local s_inv = had_pow(s, -1)
+
+        local u = sigmoid((d_sat + p_sat) * dg_sat, -1)
+        local ue = Matrix.new_vector(p_degree):fill(u)
+        local dg_aug = had(s_inv, duality_gap - ue)
+
+        local ND = diag(had(s_inv, x))
+        local N = A * ND * AT
+        local L, FD, U = M.cholesky_factorization(N)
+        L = L * FD
+
+        local r_asd = A * (-ND * dual + dg_aug) - primal
+        local y_asd = M.lu_solve_linear_equation(L, U, r_asd, fvg)
+        local s_asd = AT * -y_asd - dual
+        local x_asd = -ND * s_asd - dg_aug
+
+        -- local cor = had(x, s) + had(x, s_asd) + had(x_asd, s) + had(x_asd, s_asd)
+        -- local r_agg = Matrix.join_vector{
+        --     dual,
+        --     primal,
+        --     duality_gap + cor - cen,
+        -- }
+        -- local agg = M.gaussian_elimination(D:clone():insert_column(-r_agg), fvg(x, y, s))
+        -- local x_agg, y_agg, s_agg = split(agg)
+
+        local p_step = get_max_step(x, x_asd)
+        local d_step = get_max_step(s, s_asd)
+        debug_print(string.format(
+            "p_step = %f, d_step = %f, barrier = %f",
+            p_step, d_step, u
+        ))
+
+        x = x + p_step * x_asd
+        y = y + d_step * y_asd
+        s = s + d_step * s_asd
+    end
+    debug_print(string.format("-- complete solve %s --", problem.name))
+    debug_print("variables x:\n" .. problem:dump_primal(x))
+    debug_print("factors c:\n" .. problem:dump_primal(c))
+    -- debug_print("variables y:\n" .. problem:dump_dual(y))
+    debug_print("factors b:\n" .. problem:dump_dual(b))
+    -- debug_print("variables s:\n" .. problem:dump_primal(s))
+
+    return problem:filter_solution_to_result(x), problem:pack_pdip_variables(x, y, s)
+end
+
+--- Reduce an augmented matrix into row echelon form.
+-- @todo Refactoring for use in matrix solvers.
+-- @tparam Matrix A Matrix equation.
+-- @tparam Matrix b Column vector.
+-- @treturn Matrix Matrix of row echelon form.
+function M.gaussian_elimination(A, b)
+    local height, width = A.height, A.width
+    local ret_A = A:clone():insert_column(b)
+
+    local function select_pivot(s, x)
+        local max_value, max_index, raw_max_value = 0, nil, nil
+        for y = s, height do
+            local r = ret_A:get(y, x)
+            local a = math.abs(r)
+            if max_value < a then
+                max_value = a
+                max_index = y
+                raw_max_value = r
+            end
+        end
+        return max_index, raw_max_value
+    end
+
+    local i = 1
+    for x = 1, width + 1 do
+        local pi, pv = select_pivot(i, x)
+        if pi then
+            ret_A:row_swap(i, pi)
+            for k = i + 1, height do
+                local f = -ret_A:get(k, x) / pv
+                ret_A:row_trans(k, i, f)
+                ret_A:set(k, x, 0)
+            end
+            i = i + 1
+        end
+    end
+
+    local ret_b = ret_A:remove_column()
+    return ret_A, ret_b
+end
+
+--- LU decomposition of the symmetric matrix.
+-- @tparam Matrix A Symmetric matrix.
+-- @treturn Matrix Lower triangular matrix.
+-- @treturn Matrix Diagonal matrix.
+-- @treturn Matrix Upper triangular matrix.
+function M.cholesky_factorization(A)
+    assert(A.height == A.width)
+    local size = A.height
+    local L, D = SparseMatrix(size, size), SparseMatrix(size, size)
+    for i = 1, size do
+        local a_values = {}
+        for x, v in A:iterate_row(i) do
+            a_values[x] = v
+        end
+
+        for k = 1, i do
+            local i_it, k_it = L:iterate_row(i), L:iterate_row(k)
+            local i_r, i_v = i_it()
+            local k_r, k_v = k_it()
+
+            local sum = 0
+            while i_r and k_r do
+                if i_r < k_r then
+                    i_r, i_v = i_it()
+                elseif i_r > k_r then
+                    k_r, k_v = k_it()
+                else -- i_r == k_r
+                    local d = D:get(i_r, k_r)
+                    sum = sum + i_v * k_v * d
+                    i_r, i_v = i_it()
+                    k_r, k_v = k_it()
+                end
+            end
+
+            local a = a_values[k] or 0
+            local b = a - sum
+            if i == k then
+                D:set(k, k, math.max(b, a * machine_epsilon))
+                L:set(i, k, 1)
+            else
+                local c = D:get(k, k)
+                local v = b / c
+                L:set(i, k, v)
+            end
+        end
+    end
+    return L, D, L:T()
+end
+
+local function substitution(s, e, m, A, b, flee_value_generator)
+    local sol = {}
+    for y = s, e, m do
+        local total, factors, indexes = b:get(y, 1), {}, {}
+        for x, v in A:iterate_row(y) do
+            if sol[x] then
+                total = total - sol[x] * v
+            else
+                table.insert(factors, v)
+                table.insert(indexes, x)
+            end
+        end
+
+        local l = #indexes
+        if l == 1 then
+            sol[indexes[1]] = total / factors[1]
+        elseif l >= 2 then
+            local res = flee_value_generator(total, factors, indexes)
+            for k, x in ipairs(indexes) do
+                sol[x] = res[k]
+            end
+        end
+    end
+    return Matrix.list_to_vector(sol, A.width)
+end
+
+--- Use LU-decomposed matrices to solve linear equations.
+-- @tparam Matrix L Lower triangular matrix.
+-- @tparam Matrix U Upper triangular matrix.
+-- @tparam Matrix b Column vector.
+-- @tparam function flee_value_generator Callback function that generates the value of free variable.
+-- @treturn Matrix Solution of linear equations.
+function M.lu_solve_linear_equation(L, U, b, flee_value_generator)
+    local t = M.forward_substitution(L, b, flee_value_generator)
+    return M.backward_substitution(U, t, flee_value_generator)
+end
+
+--- Use lower triangular matrix to solve linear equations.
+-- @tparam Matrix L Lower triangular matrix.
+-- @tparam Matrix b Column vector.
+-- @tparam function flee_value_generator Callback function that generates the value of free variable.
+-- @treturn Matrix Solution of linear equations.
+function M.forward_substitution(L, b, flee_value_generator)
+    return substitution(1, L.height, 1, L, b, flee_value_generator)
+end
+
+--- Use upper triangular matrix to solve linear equations.
+-- @tparam Matrix U Upper triangular matrix.
+-- @tparam Matrix b Column vector.
+-- @tparam function flee_value_generator Callback function that generates the value of free variable.
+-- @treturn Matrix Solution of linear equations.
+function M.backward_substitution(U, b, flee_value_generator)
+    return substitution(U.height, 1, -1, U, b, flee_value_generator)
+end
+
+--- Create to callback function that generates the value of free variable.
+-- @tparam Matrix Vector to be referenced in debug output.
+-- @treturn function Callback function.
+function M.create_default_flee_value_generator(...)
+    local currents = Matrix.join_vector{...}
+    return function(target, factors, indexes)
+        debug_print(string.format("generate flee values: target = %f", target))
+        local tf = 0
+        for _, v in ipairs(factors) do
+            tf = tf + math.abs(v)
+        end
+        local ret = {}
+        local sol = target / tf
+        for i, k in ipairs(indexes) do
+            ret[i] = sol * factors[i] / math.abs(factors[i])
+            debug_print(string.format(
+                "index = %i, factor = %f, current = %f, solution = %f",
+                k, factors[i], currents[k][1], sol
+            ))
+        end
+        return ret
+    end
+end
+
+return M
\ No newline at end of file
diff --git a/modfiles/data/calculation/solver_util.lua b/modfiles/data/calculation/solver_util.lua
new file mode 100644
index 00000000..868830e9
--- /dev/null
+++ b/modfiles/data/calculation/solver_util.lua
@@ -0,0 +1,438 @@
+--- The utility for converting raw FP data into forms that can be easily processed by solvers.
+-- @module solver_util
+-- @license MIT
+-- @author B_head
+-- @todo Temperature support.
+
+local M = {}
+local calculation = calculation -- require "data.calculation.util"
+local structures = structures -- require "data.calculation.structures"
+
+local tolerance = MARGIN_OF_ERROR
+
+local function to_item_id(item_proto)
+    return item_proto.type.."@"..item_proto.name
+end
+
+local function create_item_node(item_proto, amount_per_machine_by_second)
+    local ret = {
+        normalized_id = to_item_id(item_proto),
+        name = item_proto.name,
+        type = item_proto.type,
+        amount_per_machine_by_second = amount_per_machine_by_second,
+        neighbor_recipe_lines = {}, -- Used in @{make_prioritized_links}.
+    }
+    return ret
+end
+
+local function normalize_recipe_line(recipe_line, parent)
+    local recipe_proto = recipe_line.recipe_proto
+    local machine_proto = recipe_line.machine_proto
+    local fuel_proto = recipe_line.fuel_proto
+    local total_effects = recipe_line.total_effects
+    
+    -- Not per tick.
+    local crafts_per_second = calculation.util.determine_crafts_per_tick(
+        machine_proto, recipe_proto, total_effects
+    )
+    local energy_consumption_per_machine = calculation.util.determine_energy_consumption(
+        machine_proto, 1, total_effects
+    )
+    local pollution_per_machine = calculation.util.determine_pollution(
+        machine_proto, recipe_proto, fuel_proto, total_effects, energy_consumption_per_machine
+    )
+    local production_ratio_per_machine_by_second = calculation.util.determine_production_ratio(
+        crafts_per_second, 1, 1, machine_proto.launch_sequence_time
+    )
+    local products = {}
+    for _, v in pairs(recipe_proto.products) do
+        local p_amount = calculation.util.determine_prodded_amount(
+            v, crafts_per_second, total_effects
+        )
+        local n = create_item_node(v, p_amount * crafts_per_second)
+        products[n.normalized_id] = n
+    end
+    
+    local ingredients = {}
+    for _, v in pairs(recipe_proto.ingredients) do
+        local n = create_item_node(v, v.amount * crafts_per_second)
+        ingredients[n.normalized_id] = n
+    end
+    
+    local fuel_consumption_per_machine_by_second = 0
+    if fuel_proto then
+        fuel_consumption_per_machine_by_second = calculation.util.determine_fuel_amount(
+            energy_consumption_per_machine, machine_proto.burner, fuel_proto.fuel_value, 1
+        )
+        local a = ingredients[to_item_id(fuel_proto)]
+        if a then
+            a.amount_per_machine_by_second = a.amount_per_machine_by_second + fuel_consumption_per_machine_by_second
+        else
+            local n = create_item_node(fuel_proto, fuel_consumption_per_machine_by_second)
+            ingredients[n.normalized_id] = n
+        end
+    end
+    if machine_proto.energy_type == "void" then
+        energy_consumption_per_machine = 0
+    end
+
+    local maximum_machine_count, minimum_machine_count
+    local machine_limit = recipe_line.machine_limit
+    if machine_limit then
+        maximum_machine_count = machine_limit.limit
+        if machine_limit.force_limit then
+            minimum_machine_count = machine_limit.limit
+        end
+    end
+    
+    return {
+        type = "recipe_line",
+        parent = parent,
+        floor_id = parent.floor_id,
+        line_id = recipe_line.id,
+        normalized_id = "recipe_line@" .. parent.floor_id .. "." .. recipe_line.id,
+        
+        products = products,
+        ingredients = ingredients,
+        fuel_id = fuel_proto and to_item_id(fuel_proto),
+        
+        energy_consumption_per_machine = energy_consumption_per_machine,
+        pollution_per_machine = pollution_per_machine,
+        production_ratio_per_machine_by_second = production_ratio_per_machine_by_second,
+        fuel_consumption_per_machine_by_second = fuel_consumption_per_machine_by_second,
+
+        maximum_machine_count = maximum_machine_count,
+        minimum_machine_count = minimum_machine_count,
+    }
+end
+
+local function normalize_floor(floor, line_id, parent)
+    local ret = {
+        type = "floor",
+        parent = parent,
+        floor_id = floor.id,
+        line_id = line_id,
+        normalized_id = "floor@" .. floor.id,
+        lines = {},
+    }
+    
+    for _, v in ipairs(floor.lines) do
+        local res
+        if v.subfloor then
+            res = normalize_floor(v.subfloor, v.id, ret)
+        else
+            res = normalize_recipe_line(v, ret)
+        end
+        table.insert(ret.lines, res)
+    end
+    
+    return ret
+end
+
+local function link_products_to_ingredients(from, to)
+    for _, a in pairs(from.products) do
+        for _, b in pairs(to.ingredients) do
+            if a.normalized_id == b.normalized_id then
+                table.insert(a.neighbor_recipe_lines, to)
+            end
+        end
+    end
+end
+
+local function link_ingredients_to_products(from, to)
+    for _, a in pairs(from.ingredients) do
+        for _, b in pairs(to.products) do
+            if a.normalized_id == b.normalized_id then
+                table.insert(a.neighbor_recipe_lines, to)
+            end
+        end
+    end
+end
+
+local function make_prioritized_links(normalized_top_floor)
+    for _, from in M.visit_priority_order(normalized_top_floor) do
+        for _, to in M.visit_priority_order(from.parent) do
+            link_products_to_ingredients(from, to)
+            link_ingredients_to_products(from, to)
+        end
+    end
+    return normalized_top_floor
+end
+
+--- Unify parameter units and add a true ID.
+-- Furthermore, determine the priority of item transfer.
+-- @param top_floor The value of `top_floor` in `subfactory_data`.
+-- @return Normalized data. This maintains the structure of original sub-floor.
+function M.normalize(top_floor)
+    local ret = normalize_floor(top_floor, nil, nil)
+    return make_prioritized_links(ret)
+end
+
+--- Unify parameter units of refarences.
+-- @param references The value of `top_level_products` in `subfactory_data`.
+-- @tparam number timescale The value of `timescale` in `subfactory_data`.
+-- @return Normalized references data.
+function M.normalize_references(references, timescale)
+    local ret = {}
+    for _, v in ipairs(references) do
+        local proto = v.proto
+        ret[to_item_id(proto)] = {
+            name = proto.name,
+            type = proto.type,
+            amount_per_second = v.amount / timescale
+        }
+    end
+    return ret
+end
+
+-- Variants that do not remove the element if the value is 0.
+local function class_add(class, item, amount)
+    item = (item.proto ~= nil) and item.proto or item
+    local t, n = item.type, item.name
+    class[t][n] = (class[t][n] or 0) + (amount or item.amount)
+end
+
+local function feedback_recipe_line(machine_counts, player_index, timescale, normalized_recipe_line)
+    local nrl = normalized_recipe_line
+    local machine_count = machine_counts[nrl.normalized_id]
+    
+    local energy_consumption = nrl.energy_consumption_per_machine * machine_count --- @todo Energy_consumption when idle.
+    local pollution = nrl.pollution_per_machine * machine_count --- @todo Pollution when idle.
+    local production_ratio = nrl.production_ratio_per_machine_by_second * machine_count * timescale
+    local fuel_amount = nrl.fuel_consumption_per_machine_by_second * machine_count * timescale
+    if nrl.fuel_id then
+        energy_consumption = 0
+    end
+    
+    local Product = structures.class.init()
+    for _, v in pairs(nrl.products) do
+        local amount = v.amount_per_machine_by_second * machine_count * timescale
+        class_add(Product, v, amount)
+    end
+    
+    local Ingredient = structures.class.init()
+    for k, v in pairs(nrl.ingredients) do
+        local amount = v.amount_per_machine_by_second * machine_count * timescale
+        if k == nrl.fuel_id then
+            amount = amount - fuel_amount
+            if amount > tolerance then
+                class_add(Ingredient, v, amount)
+            end
+        else
+            class_add(Ingredient, v, amount)
+        end
+    end
+    
+    return {
+        player_index = player_index,
+        floor_id = nrl.parent.floor_id,
+        line_id = nrl.line_id,
+        machine_count = machine_count,
+        energy_consumption = energy_consumption,
+        pollution = pollution,
+        production_ratio = production_ratio,
+        uncapped_production_ratio = production_ratio, --- @todo The UI code calculates and shows the change of machine count.
+        Product = Product,
+        Byproduct = structures.class.init(),
+        Ingredient = Ingredient,
+        fuel_amount = fuel_amount
+    }
+end
+
+local function class_add_all(to_class, from_class)
+    for _, item in ipairs(structures.class.to_array(from_class)) do
+        class_add(to_class, item)
+    end
+end
+
+local function class_counterbalance(class_a, class_b)
+    for _, item in ipairs(structures.class.to_array(class_b)) do
+        local t, n = item.type, item.name
+        local depot_amount = class_a[t][n] or 0
+        local counterbalance_amount = math.min(depot_amount, item.amount)
+        class_add(class_a, item, -counterbalance_amount)
+        class_add(class_b, item, -counterbalance_amount)
+    end
+end
+
+local function class_cleanup(class)
+    for _, item in ipairs(structures.class.to_array(class)) do
+        local t, n = item.type, item.name
+        if class[t][n] and class[t][n] <= tolerance then
+            class[t][n] = nil
+        end
+    end
+end
+
+local function feedback_floor(machine_counts, player_index, timescale, normalized_floor)
+    local machine_count = 0
+    local energy_consumption = 0
+    local pollution = 0
+    local Product = structures.class.init()
+    local Ingredient = structures.class.init()
+    
+    for _, l in ipairs(normalized_floor.lines) do
+        local res
+        if l.type == "floor" then
+            res = feedback_floor(machine_counts, player_index, timescale, l)
+        elseif l.type == "recipe_line" then
+            res = feedback_recipe_line(machine_counts, player_index, timescale, l)
+        else
+            assert()
+        end
+        calculation.interface.set_line_result(res)
+        
+        machine_count = machine_count + math.ceil(res.machine_count - tolerance)
+        energy_consumption = energy_consumption + res.energy_consumption
+        pollution = pollution + res.pollution
+        class_add_all(Product, res.Product)
+        class_add_all(Ingredient, res.Ingredient)
+    end
+    
+    class_counterbalance(Ingredient, Product)
+    class_cleanup(Product)
+    class_cleanup(Ingredient)
+    
+    return {
+        player_index = player_index,
+        floor_id = normalized_floor.parent and normalized_floor.parent.floor_id,
+        line_id = normalized_floor.line_id,
+        machine_count = machine_count,
+        energy_consumption = energy_consumption,
+        pollution = pollution,
+        Product = Product,
+        Byproduct = structures.class.init(),
+        Ingredient = Ingredient,
+    }
+end
+
+--- Reflect the solution in the UI.
+-- @tparam {[string]=number,...} machine_counts Solution for amount of machines needed in a sub-factory.
+-- @tparam number player_index The value of `player_index` in `subfactory_data`.
+-- @tparam number timescale The value of `timescale` in `subfactory_data`.
+-- @param normalized_top_floor The value returned by @{normalize}.
+function M.feedback(machine_counts, player_index, timescale, normalized_top_floor)
+    local res = feedback_floor(machine_counts, player_index, timescale, normalized_top_floor)
+    local ReferencesMet = structures.class.init() --- @todo Reflect on results.
+    calculation.interface.set_subfactory_result{
+        player_index = player_index,
+        energy_consumption = res.energy_consumption,
+        pollution = res.pollution,
+        Product = ReferencesMet,
+        Byproduct = res.Product,
+        Ingredient = res.Ingredient
+    }
+end
+
+-- No coroutine? Damn it!
+--
+-- local function icoroutine(func, a1, a2, a3, a4, a5, a6, a7 ,a8, a9)
+--     local success
+--     return function(co)
+--         success, a1, a2, a3, a4, a5, a6, a7 ,a8, a9 = coroutine.resume(co, a1, a2, a3, a4, a5, a6, a7 ,a8, a9)
+--         if not success then
+--             error(debug.traceback(co, a1), 2)
+--         end
+--         return a1, a2, a3, a4, a5, a6, a7 ,a8, a9
+--     end, coroutine.create(func)
+-- end
+--
+-- local function visit_priority_order(floor, prev_floor)
+--     for _, v in ipairs(floor.lines) do
+--         if v.type == "floor" then
+--             if v ~= prev_floor then
+--                 visit_priority_order(v, floor)
+--             end
+--         elseif v.type == "recipe_line" then
+--             coroutine.yield(v)
+--         else
+--             assert(false)
+--         end
+--     end
+--     local parent = floor.parent
+--     if parent and parent ~= prev_floor then
+--         visit_priority_order(parent, floor)
+--     end
+-- end
+
+--- Visit the recipe lines in the following order.
+--
+-- 1. Visit the recipe lines in start_floor from top to bottom.
+-- 2. Visit the sub-floor from top to bottom, recursively, depth-first.
+--    Within each sub-floor, visit the recipe lines from top to bottom.
+-- 3. Visit the super-floor, then visit the recipe lines and sub-floors in the same order as above.
+--    After that, recursively visit even higher level super-floor.
+--
+-- (I can't explain this well, sorry.)
+--
+-- @param start_floor Floor to start the visit.
+-- @return Iterator object that returns recipe lines.
+function M.visit_priority_order(start_floor)
+    local floor_stack = {start_floor}
+    local index_stack = {1}
+    local mode_stack = {"recipe_line"}
+    local function it()
+        local top = #floor_stack
+        if top == 0 then
+            return nil
+        end
+        local floor, prev_floor = floor_stack[top], floor_stack[top - 1]
+        local index, mode = index_stack[top], mode_stack[top]
+        index_stack[top] = index + 1
+        if index <= #floor.lines then
+            local v = floor.lines[index]
+            if v.type == "floor" and mode == "floor" then
+                if v ~= prev_floor then
+                    table.insert(floor_stack, v)
+                    table.insert(index_stack, 1)
+                    table.insert(mode_stack, "recipe_line")
+                end
+            elseif v.type == "recipe_line" and mode == "recipe_line"  then
+                return v.normalized_id, v
+            end
+        else
+            if mode == "recipe_line" then
+                mode_stack[top] = "floor"
+                index_stack[top] = 1
+            elseif mode == "floor" then
+                mode_stack[top] = "parent"
+                local parent = floor.parent
+                if parent and parent ~= prev_floor then
+                    table.insert(floor_stack, parent)
+                    table.insert(index_stack, 1)
+                    table.insert(mode_stack, "recipe_line")
+                end
+            elseif mode == "parent" then
+                table.remove(floor_stack)
+                table.remove(index_stack)
+                table.remove(mode_stack)
+            end
+        end
+        return it()
+    end
+    return it
+end
+
+--- Returns an iterator depending on the container structure.
+-- @param recipe_lines Container for recipe line.
+-- @return Iterator object that returns recipe lines.
+function M.iterate_recipe_lines(recipe_lines)
+    if recipe_lines.type == "floor" then
+        return M.visit_priority_order(recipe_lines)
+    else
+        return pairs(recipe_lines) -- by flat_recipe_lines.
+    end
+end
+
+--- Convert subfactories to list format.
+-- @param normalized_top_floor The value returned by @{normalize}.
+-- @return List of recipe lines.
+function M.to_flat_recipe_lines(normalized_top_floor)
+    local ret = {}
+    for k, v in M.visit_priority_order(normalized_top_floor) do
+        ret[k] = v
+    end
+    return ret
+end
+
+return M
diff --git a/modfiles/data/classes/Subfactory.lua b/modfiles/data/classes/Subfactory.lua
index fc4925a5..1637981f 100644
--- a/modfiles/data/classes/Subfactory.lua
+++ b/modfiles/data/classes/Subfactory.lua
@@ -14,7 +14,9 @@ function Subfactory.init(name, icon)
         Byproduct = Collection.init("Item"),
         Ingredient = Collection.init("Item"),
         Floor = Collection.init("Floor"),
+        solver_type = nil,  -- needs to be set after init
         matrix_free_items = nil,
+        prev_raw_solution = nil,
         linearly_dependant = false,  -- determined by the solver
         selected_floor = nil,
         item_request_proxy = nil,
@@ -169,6 +171,7 @@ function Subfactory.pack(self)
         notes = self.notes,
         mining_productivity = self.mining_productivity,
         Product = Collection.pack(self.Product),
+        solver_type = self.solver_type,
         matrix_free_items = packed_free_items,
         -- Floors get packed by recursive nesting, which is necessary for a json-type data
         -- structure. It will need to be unpacked into the regular structure 'manually'.
@@ -184,6 +187,7 @@ function Subfactory.unpack(packed_self)
     self.notes = packed_self.notes
     self.mining_productivity = packed_self.mining_productivity
     self.Product = Collection.unpack(packed_self.Product, self)
+    self.solver_type = packed_self.solver_type
 
     if packed_self.matrix_free_items then
         self.matrix_free_items = {}
diff --git a/modfiles/data/handlers/migrator.lua b/modfiles/data/handlers/migrator.lua
index b105ab09..501d5d0f 100644
--- a/modfiles/data/handlers/migrator.lua
+++ b/modfiles/data/handlers/migrator.lua
@@ -27,6 +27,7 @@ local migration_masterlist = {
     [17] = {version="1.1.25", migration=require("data.migrations.migration_1_1_25")},
     [18] = {version="1.1.26", migration=require("data.migrations.migration_1_1_26")},
     [19] = {version="1.1.27", migration=require("data.migrations.migration_1_1_27")},
+    [20] = {version="1.1.28", migration=require("data.migrations.migration_1_1_28")},
 }
 
 -- ** LOCAL UTIL **
diff --git a/modfiles/data/init.lua b/modfiles/data/init.lua
index f1f68d32..6b1eb11f 100755
--- a/modfiles/data/init.lua
+++ b/modfiles/data/init.lua
@@ -35,7 +35,7 @@ local function reload_settings(player)
     settings_table.alt_action = settings["fp_alt_action"].value
     settings_table.default_timescale = timescale_to_number[settings["fp_default_timescale"].value]
     settings_table.belts_or_lanes = settings["fp_view_belts_or_lanes"].value
-    settings_table.prefer_matrix_solver = settings["fp_prefer_matrix_solver"].value
+    settings_table.default_solver_type = settings["fp_default_solver_type"].value
 
     global.players[player.index].settings = settings_table
 end
diff --git a/modfiles/data/migrations/masterlist.json b/modfiles/data/migrations/masterlist.json
index faea1dec..1b084285 100644
--- a/modfiles/data/migrations/masterlist.json
+++ b/modfiles/data/migrations/masterlist.json
@@ -17,5 +17,6 @@
     "1.1.21",
     "1.1.25",
     "1.1.26",
-    "1.1.27"
+    "1.1.27",
+    "1.1.28"
 ]
diff --git a/modfiles/data/migrations/migration_1_1_28.lua b/modfiles/data/migrations/migration_1_1_28.lua
new file mode 100644
index 00000000..d4ed8c47
--- /dev/null
+++ b/modfiles/data/migrations/migration_1_1_28.lua
@@ -0,0 +1,25 @@
+local migration = {}
+
+function migration.global()
+end
+
+function migration.player_table(player_table)
+end
+
+function migration.subfactory(subfactory)
+    if subfactory.matrix_free_items then
+        subfactory.solver_type = "matrix"
+    else
+        subfactory.solver_type = "traditional"
+    end
+end
+
+function migration.packed_subfactory(packed_subfactory)
+    if packed_subfactory.matrix_free_items then
+        packed_subfactory.solver_type = "matrix"
+    else
+        packed_subfactory.solver_type = "traditional"
+    end
+end
+
+return migration
diff --git a/modfiles/locale/en/config.cfg b/modfiles/locale/en/config.cfg
index d2d2fc62..22b15b40 100644
--- a/modfiles/locale/en/config.cfg
+++ b/modfiles/locale/en/config.cfg
@@ -31,7 +31,7 @@ fp_subfactory_list_rows=Interface height [img=info]
 fp_alt_action=Open in other mod [img=info]
 fp_default_timescale=Default timescale [img=info]
 fp_view_belts_or_lanes=Belts or Lanes [img=info]
-fp_prefer_matrix_solver=Prefer matrix solver [img=info]
+fp_default_solver_type=Default solver [img=info]
 
 [mod-setting-description]
 fp_display_gui_button=Shows the button on the top left of the screen. It opens and closes the main interface.
@@ -40,7 +40,7 @@ fp_subfactory_list_rows=Set the main interface height by choosing how many subfa
 fp_alt_action=Select an action when alt-clicking any item- or recipe-button. Currently works with FNEI, WIIRUF and Recipe Book. The mod has to be installed and enabled before it shows up in this list.
 fp_default_timescale=Choose the timescale that any new subfactory should be created with.
 fp_view_belts_or_lanes=Indicate whether you think of item throughput as individual lanes or as full belts.
-fp_prefer_matrix_solver=Decide whether new subfactories should enable the matrix solver instead of the traditional one.
+fp_default_solver_type=Choose the solver that any new subfactory should be created with.
 
 [string-mod-setting]
 fp_alt_action-none=None
@@ -52,6 +52,9 @@ fp_default_timescale-one_minute=1 Minute
 fp_default_timescale-one_hour=1 Hour
 fp_view_belts_or_lanes-belts=Belts
 fp_view_belts_or_lanes-lanes=Lanes
+fp_default_solver_type-traditional=Traditional
+fp_default_solver_type-matrix=Matrix
+fp_default_solver_type-interior_point=Interior point
 
 
 [shortcut-name]
@@ -296,9 +299,10 @@ mining_productivity=Mining productivity
 mining_productivity_tt=The current mining productivity bonus. By default, it tracks your research progress, but you can override the percentage manually. To revert to automated tracking, clear the textfield and confirm.
 override=Override
 solver_choice=Solver
-solver_choice_tt=Choose which of the solvers to use for this subfactory. The traditional one works by going through your recipes in order and figuring out their needs. The matrix solver on the other hand can deal with loops and byproducts, but sometimes needs additional configuration.
+solver_choice_tt=Choose which of the solvers to use for this subfactory.\nThe traditional solver will visit the list in order from top, and solve to satisfy the previous needs.\nThe matrix solver supports loops, and using byproducts, but sometimes needs additional configuration.\nThe interior point solver supports loops, and allows bottom-up design, and respects the list in order and sub-floor structure like the traditional solver.
 solver_choice_traditional=Traditional
 solver_choice_matrix=Matrix
+solver_choice_interior_point=Interior point
 solver_choice_configure=Configure the matrix solver
 subfactory_modset_changes=Your active mods changed:
 subfactory_mod_removed=\n\n[color=#FF3333][font=default-bold]These mods were removed:[/font][/color]
@@ -358,6 +362,7 @@ tut_mode_machine_matrix=Left-click: Change machine\nShift-left-click: Upgrade\nC
 tut_mode_beacon=Left/Right-click: Edit\nControl-right-click: Delete\nAlt-left-click: Put into cursor
 tut_mode_module=Left/Right-click: Edit\nControl-right-click: Delete
 tut_mode_product=Left-click: Prioritize\nRight-click: Specify amount
+tut_mode_product_matrix=Left-click: Add recipe to the end\nShift-left-click: Add recipe right below
 tut_mode_byproduct=Right-click: Specify amount
 tut_mode_byproduct_matrix=Left-click: Add recipe to the end\nShift-left-click: Add recipe right below
 tut_mode_ingredient=Left-click: Add recipe to the end\nRight-click: Specify amount\nShift-left-click: Add recipe right below
@@ -493,6 +498,8 @@ pl_floor=__plural_for_parameter_1_{1=floor|rest=floors}__
 pu_floor=__plural_for_parameter_1_{1=Floor|rest=Floors}__
 pl_item=__plural_for_parameter_1_{1=item|rest=items}__
 pu_item=__plural_for_parameter_1_{1=Item|rest=Items}__
+pl_refarence=__plural_for_parameter_1_{1=refarence|rest=refarences}__
+pu_refarence=__plural_for_parameter_1_{1=Refarence|rest=Refarences}__
 pl_product=__plural_for_parameter_1_{1=product|rest=products}__
 pu_product=__plural_for_parameter_1_{1=Product|rest=Products}__
 pl_byproduct=__plural_for_parameter_1_{1=byproduct|rest=byproducts}__
diff --git a/modfiles/settings.lua b/modfiles/settings.lua
index c1ece8a8..8461d875 100644
--- a/modfiles/settings.lua
+++ b/modfiles/settings.lua
@@ -59,10 +59,11 @@ data:extend({
         order = "f"
     },
     {
-        type = "bool-setting",
-        name = "fp_prefer_matrix_solver",
+        type = "string-setting",
+        name = "fp_default_solver_type",
         setting_type = "runtime-per-user",
-        default_value = false,
+        default_value = "interior_point",
+        allowed_values = {"traditional", "matrix", "interior_point"},
         order = "g"
-    },
+    }
 })
diff --git a/modfiles/ui/elements/item_boxes.lua b/modfiles/ui/elements/item_boxes.lua
index 0c2fd51c..c6cf49b9 100644
--- a/modfiles/ui/elements/item_boxes.lua
+++ b/modfiles/ui/elements/item_boxes.lua
@@ -2,7 +2,7 @@ item_boxes = {}
 
 --- ** LOCAL UTIL **
 local function add_recipe(player, context, type, item_proto)
-    if type == "byproduct" and context.subfactory.matrix_free_items == nil then
+    if type == "byproduct" and context.subfactory.solver_type == "traditional" then
         title_bar.enqueue_message(player, {"fp.error_cant_add_byproduct_recipe"}, "error", 1, true)
         return
     end
@@ -21,14 +21,15 @@ end
 local function build_item_box(player, category, column_count)
     local item_boxes_elements = data_util.get("main_elements", player).item_boxes
 
-    local window_frame = item_boxes_elements.horizontal_flow.add{type="frame", direction="vertical",
-      style="inside_shallow_frame"}
+    local window_frame = item_boxes_elements.horizontal_flow.add{type="frame",
+      direction="vertical", style="inside_shallow_frame"}
     window_frame.style.top_padding = 6
     window_frame.style.bottom_padding = ITEM_BOX_PADDING
 
-    local label = window_frame.add{type="label", caption={"fp.pu_" .. category, 2}, style="caption_label"}
+    local label = window_frame.add{type="label", style="caption_label"}
     label.style.left_padding = ITEM_BOX_PADDING
     label.style.bottom_margin = 4
+    item_boxes_elements[category .. "_category_name"] = label
 
     local scroll_pane = window_frame.add{type="scroll-pane", style="fp_scroll-pane_slot_table"}
     scroll_pane.style.maximal_height = ITEM_BOX_MAX_ROWS * ITEM_BOX_BUTTON_SIZE
@@ -55,13 +56,22 @@ local function refresh_item_box(player, category, subfactory, allow_addition)
     local table_item_count = 0
     local metadata = view_state.generate_metadata(player, subfactory, 4, true)
     local default_style, tut_mode_tooltip, global_enable = "flib_slot_button_default", "", true
+    local solver_type = subfactory.solver_type
+    
+    local label = item_boxes_elements[category .. "_category_name"]
+    if subfactory.solver_type == "interior_point" then
+        local conv = {["product"]="refarence", ["byproduct"]="product", ["ingredient"]="ingredient"}
+        label.caption = {"fp.pu_" .. conv[category], 2}
+    else
+        label.caption = {"fp.pu_" .. category, 2}
+    end
 
     if class == "Product" then
         default_style = "flib_slot_button_red"
         tut_mode_tooltip = ui_util.generate_tutorial_tooltip(player, "tl_product", true, true, true)
     elseif class == "Byproduct" then
         default_style = "flib_slot_button_red"
-        if subfactory.matrix_free_items ~= nil then
+        if solver_type == "matrix" or solver_type == "interior_point" then
             tut_mode_tooltip = ui_util.generate_tutorial_tooltip(player, "tl_byproduct", true, true, true)
         else
             global_enable = false
diff --git a/modfiles/ui/elements/production_handler.lua b/modfiles/ui/elements/production_handler.lua
index 1c195bfc..f46ad6bb 100644
--- a/modfiles/ui/elements/production_handler.lua
+++ b/modfiles/ui/elements/production_handler.lua
@@ -200,6 +200,7 @@ local function handle_machine_click(player, tags, metadata)
     if not ui_util.check_archive_status(player) then return end
 
     local context = data_util.get("context", player)
+    local solver_type = context.subfactory.solver_type
     local line = Floor.get(context.floor, "Line", tags.line_id)
     -- I don't need to care about relevant lines here because this only gets called on lines without subfloor
 
@@ -265,7 +266,7 @@ local function handle_machine_click(player, tags, metadata)
             calculation.update(player, context.subfactory)
             main_dialog.refresh(player, "subfactory")
 
-        elseif context.subfactory.matrix_free_items == nil then  -- open the machine options
+        elseif solver_type == "traditional" or solver_type == "interior_point" then  -- open the machine limit options
             local modal_data = {
                 title = {"fp.options_machine_title"},
                 text = {"fp.options_machine_text", line.machine.proto.localised_name},
@@ -378,6 +379,7 @@ end
 
 local function handle_item_click(player, tags, metadata)
     local context = data_util.get("context", player)
+    local solver_type = context.subfactory.solver_type
     local line = Floor.get(context.floor, "Line", tags.line_id)
     -- I don't need to care about relevant lines here because this only gets called on lines without subfloor
     local item = Line.get(line, tags.class, tags.item_id)
@@ -390,52 +392,64 @@ local function handle_item_click(player, tags, metadata)
 
     elseif metadata.click == "left" and item.proto.type ~= "entity" then  -- Handles the specific type of item actions
         if tags.class == "Product" then -- Set the priority product
-            if line.Product.count < 2 then
-                title_bar.enqueue_message(player, {"fp.warning_no_prioritizing_single_product"}, "warning", 1, true)
-            elseif context.subfactory.matrix_free_items == nil then
-                -- Remove the priority_product if the already selected one is clicked
-                line.priority_product_proto = (line.priority_product_proto ~= item.proto) and item.proto or nil
-
-                calculation.update(player, context.subfactory)
-                main_dialog.refresh(player, "subfactory")
+            if solver_type == "traditional" then
+                if line.Product.count < 2 then
+                    title_bar.enqueue_message(player, {"fp.warning_no_prioritizing_single_product"}, "warning", 1, true)
+                else
+                    -- Remove the priority_product if the already selected one is clicked
+                    line.priority_product_proto = (line.priority_product_proto ~= item.proto) and item.proto or nil
+
+                    calculation.update(player, context.subfactory)
+                    main_dialog.refresh(player, "subfactory")
+                end
+            elseif solver_type == "interior_point" then
+                modal_dialog.enter(player, {type="recipe", modal_data={product_proto=item.proto,
+                  production_type="consume", add_after_position=((metadata.shift) and line.gui_position or nil)}})
             end
-
-        else  -- Byproduct or Ingredient
-            local production_type = (tags.class == "Byproduct") and "consume" or "produce"
-            -- The sequential solver does not support byproduct recipes at the moment
-            if production_type == "consume" and context.subfactory.matrix_free_items == nil then
+        elseif tags.class == "Byproduct" then
+            if solver_type == "traditional" then
                 title_bar.enqueue_message(player, {"fp.error_cant_add_byproduct_recipe"}, "error", 1, true)
             else
                 modal_dialog.enter(player, {type="recipe", modal_data={product_proto=item.proto,
-                  production_type=production_type, add_after_position=((metadata.shift) and line.gui_position or nil)}})
+                  production_type="consume", add_after_position=((metadata.shift) and line.gui_position or nil)}})
             end
+        elseif tags.class == "Ingredient" then
+            modal_dialog.enter(player, {type="recipe", modal_data={product_proto=item.proto,
+                production_type="produce", add_after_position=((metadata.shift) and line.gui_position or nil)}})
+        else
+            assert()
         end
 
-    elseif metadata.click == "right" and context.subfactory.matrix_free_items == nil then
-        -- Set the view state so that the amount shown in the dialog makes sense
-        view_state.select(player, "items_per_timescale", "subfactory")  -- refreshes "subfactory" if necessary
+    elseif metadata.click == "right" then
+        if solver_type == "traditional" then
+            -- Set the view state so that the amount shown in the dialog makes sense
+            view_state.select(player, "items_per_timescale", "subfactory")  -- refreshes "subfactory" if necessary
 
-        local type_localised_string = {"fp.pl_" .. tags.class:lower(), 1}
-        local produce_consume = (tags.class == "Ingredient") and {"fp.consume"} or {"fp.produce"}
+            local type_localised_string = {"fp.pl_" .. tags.class:lower(), 1}
+            local produce_consume = (tags.class == "Ingredient") and {"fp.consume"} or {"fp.produce"}
 
-        local modal_data = {
-            title = {"fp.options_item_title", type_localised_string},
-            text = {"fp.options_item_text", item.proto.localised_name},
-            submission_handler_name = "apply_item_options",
-            object = item,
-            fields = {
-                {
-                    type = "numeric_textfield",
-                    name = "item_amount",
-                    caption = {"fp.options_item_amount"},
-                    tooltip = {"fp.options_item_amount_tt", type_localised_string, produce_consume},
-                    text = item.amount,
-                    width = 140,
-                    focus = true
+            local modal_data = {
+                title = {"fp.options_item_title", type_localised_string},
+                text = {"fp.options_item_text", item.proto.localised_name},
+                submission_handler_name = "apply_item_options",
+                object = item,
+                fields = {
+                    {
+                        type = "numeric_textfield",
+                        name = "item_amount",
+                        caption = {"fp.options_item_amount"},
+                        tooltip = {"fp.options_item_amount_tt", type_localised_string, produce_consume},
+                        text = item.amount,
+                        width = 140,
+                        focus = true
+                    }
                 }
             }
-        }
-        modal_dialog.enter(player, {type="options", modal_data=modal_data})
+            modal_dialog.enter(player, {type="options", modal_data=modal_data})
+        elseif solver_type == "interior_point" then
+            -- todo: Open the machine limit and item limit options in the improved dialog.
+        end
+        
     end
 end
 
diff --git a/modfiles/ui/elements/production_table.lua b/modfiles/ui/elements/production_table.lua
index 635daa9a..9469545e 100644
--- a/modfiles/ui/elements/production_table.lua
+++ b/modfiles/ui/elements/production_table.lua
@@ -11,7 +11,7 @@ local function generate_metadata(player)
 
     local metadata = {
         archive_open = (ui_state.flags.archive_open),
-        matrix_solver_active = (subfactory.matrix_free_items ~= nil),
+        solver_type = subfactory.solver_type,
         mining_productivity = mining_productivity,
         round_button_numbers = preferences.round_button_numbers,
         pollution_column = preferences.pollution_column,
@@ -22,7 +22,32 @@ local function generate_metadata(player)
 
     if preferences.tutorial_mode then
         -- Choose the right type of tutorial text right here if possible
-        local matrix_postfix = (metadata.matrix_solver_active) and "_matrix" or ""
+        local element_types
+        if metadata.solver_type == "traditional" then
+          element_types = {
+            ["machine"] = "machine",
+            ["product"] = "product",
+            ["byproduct"] = "byproduct",
+            ["ingredient"] = "ingredient",
+            ["ingredient_entity"] = "ingredient_entity",
+          }
+        elseif metadata.solver_type == "matrix" then
+          element_types = {
+            ["machine"] = "machine_matrix",
+            ["product"] = nil,
+            ["byproduct"] = "byproduct_matrix",
+            ["ingredient"] = "ingredient_matrix",
+            ["ingredient_entity"] = nil,
+          }
+        elseif metadata.solver_type == "interior_point" then
+          element_types = {
+            ["machine"] = "machine",
+            ["product"] = "product_matrix",
+            ["byproduct"] = nil,
+            ["ingredient"] = "ingredient_matrix",
+            ["ingredient_entity"] = nil,
+          }
+        end
 
         metadata.production_toggle_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "production_toggle",
           false, false, true)
@@ -30,16 +55,16 @@ local function generate_metadata(player)
           true, true, true)
         metadata.consuming_recipe_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "consuming_recipe",
           true, true, true)
-        metadata.machine_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "machine" .. matrix_postfix,
+        metadata.machine_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, element_types["machine"],
           false, true, true)
         metadata.beacon_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "beacon", false, true, true)
         metadata.module_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "module", false, true, true)
-        metadata.product_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "product", true, true, true)
-        metadata.byproduct_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "byproduct" .. matrix_postfix,
+        metadata.product_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, element_types["product"], true, true, true)
+        metadata.byproduct_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, element_types["byproduct"],
          true, true, true)
-        metadata.ingredient_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "ingredient" .. matrix_postfix,
+        metadata.ingredient_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, element_types["ingredient"],
           true, true, true)
-        metadata.ingredient_entity_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "ingredient_entity",
+        metadata.ingredient_entity_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, element_types["ingredient_entity"],
           true, true, true)
         metadata.fuel_tutorial_tooltip = ui_util.generate_tutorial_tooltip(player, "fuel", true, true, true)
     end
@@ -97,7 +122,7 @@ end
 function builders.percentage(line, parent_flow, metadata)
     local relevant_line = (line.subfloor) and line.subfloor.defining_line or line
 
-    local enabled = (not metadata.archive_open and not metadata.matrix_solver_active)
+    local enabled = (not metadata.archive_open)
     local textfield_percentage = parent_flow.add{type="textfield", text=tostring(relevant_line.percentage),
       tags={mod="fp", on_gui_text_changed="line_percentage", on_gui_confirmed="line_percentage", line_id=line.id},
       enabled=enabled}
@@ -126,7 +151,7 @@ function builders.machine(line, parent_flow, metadata)
         if metadata.round_button_numbers then machine_count = math.ceil(machine_count) end
 
         local style, indication, machine_limit = "flib_slot_button_default_small", "", line.machine.limit
-        if not metadata.matrix_solver_active and machine_limit ~= nil then
+        if metadata.solver_type ~= "matrix" and machine_limit ~= nil then
             if line.machine.force_limit then
                 style = "flib_slot_button_pink_small"
                 indication = {"fp.newline", {"fp.notice", {"fp.machine_limit_force", machine_limit}}}
@@ -246,18 +271,22 @@ function builders.products(line, parent_flow, metadata)
         -- items/s/machine does not make sense for lines with subfloors, show items/s instead
         local machine_count = (not line.subfloor) and line.machine.count or nil
         local amount, number_tooltip = view_state.process_item(metadata.view_state_metadata,
-          product, nil, machine_count)
+          product, nil, machine_count, not line.subfloor)
         if amount == "0" and line.subfloor then goto skip_product end  -- amount can't be -1 for products
 
-        local style = "flib_slot_button_default_small"
+        local style, enabled = "flib_slot_button_default_small", false
         local indication_string, tutorial_tooltip = "", ""
 
-        if not line.subfloor and not metadata.matrix_solver_active then
+        if metadata.solver_type == "traditional" and not line.subfloor then
             -- We can check for identity because they reference the same table
             if line.priority_product_proto == product.proto then
                 style = "flib_slot_button_pink_small"
                 indication_string = {"fp.indication", {"fp.priority_product"}}
             end
+            enabled = true
+            tutorial_tooltip = metadata.product_tutorial_tooltip
+        else
+            enabled = (metadata.solver_type ~= "matrix")
             tutorial_tooltip = metadata.product_tutorial_tooltip
         end
 
@@ -267,7 +296,7 @@ function builders.products(line, parent_flow, metadata)
 
         parent_flow.add{type="sprite-button", tags={mod="fp", on_gui_click="act_on_line_item", line_id=line.id,
           class="Product", item_id=product.id}, sprite=product.proto.sprite, style=style, number=amount,
-          tooltip=tooltip, enabled=(not line.subfloor), mouse_button_filter={"left-and-right"}}
+          tooltip=tooltip, enabled=enabled, mouse_button_filter={"left-and-right"}}
 
         ::skip_product::
     end
@@ -278,11 +307,12 @@ function builders.byproducts(line, parent_flow, metadata)
         -- items/s/machine does not make sense for lines with subfloors, show items/s instead
         local machine_count = (not line.subfloor) and line.machine.count or nil
         local amount, number_tooltip = view_state.process_item(metadata.view_state_metadata,
-          byproduct, nil, machine_count)
+          byproduct, nil, machine_count, not line.subfloor)
         if amount == -1 then goto skip_byproduct end  -- an amount of -1 means it was below the margin of error
 
+        local enabled = (metadata.solver_type ~= "traditional" or not line.subfloor)
         local number_line = (number_tooltip) and {"fp.newline", number_tooltip} or ""
-        local tutorial_tooltip = (not line.subfloor) and metadata.byproduct_tutorial_tooltip or ""
+        local tutorial_tooltip = enabled and metadata.byproduct_tutorial_tooltip or ""
         local tooltip = {"", byproduct.proto.localised_name, number_line, tutorial_tooltip}
 
         parent_flow.add{type="sprite-button", tags={mod="fp", on_gui_click="act_on_line_item", line_id=line.id,
@@ -298,7 +328,7 @@ function builders.ingredients(line, parent_flow, metadata)
         -- items/s/machine does not make sense for lines with subfloors, show items/s instead
         local machine_count = (not line.subfloor) and line.machine.count or nil
         local amount, number_tooltip = view_state.process_item(metadata.view_state_metadata,
-          ingredient, nil, machine_count)
+          ingredient, nil, machine_count, not line.subfloor)
         if amount == -1 then goto skip_ingredient end  -- an amount of -1 means it was below the margin of error
 
         local style, enabled = "flib_slot_button_green_small", true
@@ -307,11 +337,11 @@ function builders.ingredients(line, parent_flow, metadata)
 
         if ingredient.proto.type == "entity" then
             style = "flib_slot_button_default_small"
-            enabled = (not metadata.matrix_solver_active)
+            enabled = (metadata.solver_type == "traditional")
             indication_string = {"fp.indication", {"fp.raw_ore"}}
-            tutorial_tooltip = (metadata.matrix_solver_active) and "" or metadata.ingredient_entity_tutorial_tooltip
+            tutorial_tooltip = metadata.ingredient_entity_tutorial_tooltip
 
-        elseif metadata.ingredient_satisfaction then
+        elseif metadata.ingredient_satisfaction and ingredient.amount > 0 then
             local satisfaction_percentage = (ingredient.satisfied_amount / ingredient.amount) * 100
             local formatted_percentage = ui_util.format_number(satisfaction_percentage, 3)
 
@@ -344,7 +374,8 @@ end
 function builders.fuel(line, parent_flow, metadata)
     local fuel = line.machine.fuel
 
-    local amount, number_tooltip = view_state.process_item(metadata.view_state_metadata, fuel, nil, line.machine.count)
+    local amount, number_tooltip = view_state.process_item(metadata.view_state_metadata,
+      fuel, nil, line.machine.count, not line.subfloor)
     if amount == -1 then return end  -- an amount of -1 means it was below the margin of error
 
     local satisfaction_line = ""
@@ -414,10 +445,18 @@ function production_table.refresh(player)
     production_table_elements.production_scroll_pane.visible = (subfactory_valid and any_lines_present)
     if not subfactory_valid then return end
 
+    local show_column = {}
+    if subfactory.solver_type == "interior_point" then
+        show_column = {["percentage"]=false, ["byproducts"]=false}
+    elseif subfactory.solver_type == "matrix" then
+        show_column = {["percentage"]=false}
+    end
+
     local production_columns, column_count = {}, 0
     for _, column_data in ipairs(all_production_columns) do
         -- Explicit comparison needed here, as both true and nil columns should be shown
-        if preferences[column_data.name .. "_column"] ~= false then
+        local name = column_data.name
+        if show_column[name] ~= false and preferences[name .. "_column"] ~= false then
             column_count = column_count + 1
             production_columns[column_count] = column_data
         end
diff --git a/modfiles/ui/elements/subfactory_info.lua b/modfiles/ui/elements/subfactory_info.lua
index a29a0532..6480e99b 100644
--- a/modfiles/ui/elements/subfactory_info.lua
+++ b/modfiles/ui/elements/subfactory_info.lua
@@ -35,13 +35,18 @@ end
 
 local function handle_solver_change(player, _, metadata)
     local subfactory = data_util.get("context", player).subfactory
-    local new_solver = (metadata.switch_state == "left") and "traditional" or "matrix"
+    local new_solver = SOLVER_TYPE_MAP[metadata.selected_index]
+
+    subfactory.solver_type = new_solver
+    subfactory.prev_raw_solution = nil
+    subfactory.linearly_dependant = false
 
     if new_solver == "matrix" then
-        subfactory.matrix_free_items = {}  -- 'activate' the matrix solver
-    else
-        subfactory.matrix_free_items = nil  -- disable the matrix solver
-        subfactory.linearly_dependant = false
+        subfactory.matrix_free_items = {}
+    elseif new_solver == "interior_point" then
+        subfactory.matrix_free_items = nil
+    elseif new_solver == "traditional" then
+        subfactory.matrix_free_items = nil
 
         -- This function works its way through subfloors. Consuming recipes can't have subfloors though.
         local any_lines_removed = false
@@ -63,6 +68,8 @@ local function handle_solver_change(player, _, metadata)
         if any_lines_removed then  -- inform the user if any byproduct recipes are being removed
             title_bar.enqueue_message(player, {"fp.hint_byproducts_removed"}, "hint", 1, false)
         end
+    else
+        assert(false, new_solver)
     end
 
     calculation.update(player, subfactory)
@@ -195,10 +202,14 @@ function subfactory_info.build(player)
     flow_solver_choice.add{type="label", caption={"fp.key_title", {"fp.info_label", {"fp.solver_choice"}}},
       tooltip={"fp.solver_choice_tt"}}
 
-    local switch_solver_choice = flow_solver_choice.add{type="switch", right_label_caption={"fp.solver_choice_matrix"},
-      left_label_caption={"fp.solver_choice_traditional"},
-      tags={mod="fp", on_gui_switch_state_changed="solver_choice_changed"}}
-    main_elements.subfactory_info["solver_choice_switch"] = switch_solver_choice
+    local dropdown_solver_choice = flow_solver_choice.add{type="drop-down",
+      items={
+        {"fp.solver_choice_traditional"},
+        {"fp.solver_choice_matrix"},
+        {"fp.solver_choice_interior_point"},
+      },
+      tags={mod="fp", on_gui_selection_state_changed="solver_choice_changed"}}
+    main_elements.subfactory_info["solver_choice_dropdown"] = dropdown_solver_choice
 
     local button_configure_solver = flow_solver_choice.add{type="sprite-button", sprite="utility/change_recipe",
       tooltip={"fp.solver_choice_configure"}, tags={mod="fp", on_gui_click="configure_matrix_solver"},
@@ -231,7 +242,6 @@ function subfactory_info.refresh(player)
 
     elseif valid_subfactory_selected then  -- we need to refresh some stuff in this case
         local archive_open = ui_state.flags.archive_open
-        local matrix_solver_active = (subfactory.matrix_free_items ~= nil)
 
         -- Power + Pollution
         local label_power = subfactory_info_elements.power_label
@@ -281,10 +291,14 @@ function subfactory_info.refresh(player)
         subfactory_info_elements.percentage_label.visible = custom_prod_set
 
         -- Solver Choice
-        local switch_state = (matrix_solver_active) and "right" or "left"
-        subfactory_info_elements.solver_choice_switch.switch_state = switch_state
-        subfactory_info_elements.solver_choice_switch.enabled = (not archive_open)
-        subfactory_info_elements.configure_solver_button.enabled = (not archive_open and matrix_solver_active)
+        local solver_type = subfactory.solver_type
+        for index, value in ipairs(SOLVER_TYPE_MAP) do
+            if solver_type == value then
+                subfactory_info_elements.solver_choice_dropdown.selected_index = index
+            end
+        end
+        subfactory_info_elements.solver_choice_dropdown.enabled = (not archive_open)
+        subfactory_info_elements.configure_solver_button.enabled = (not archive_open and solver_type == "matrix")
     end
 end
 
@@ -335,7 +349,7 @@ subfactory_info.gui_events = {
             end)
         }
     },
-    on_gui_switch_state_changed = {
+    on_gui_selection_state_changed = {
         {
             name = "solver_choice_changed",
             handler = handle_solver_change
diff --git a/modfiles/ui/elements/subfactory_list.lua b/modfiles/ui/elements/subfactory_list.lua
index 775019e9..cdda2cc6 100644
--- a/modfiles/ui/elements/subfactory_list.lua
+++ b/modfiles/ui/elements/subfactory_list.lua
@@ -339,7 +339,7 @@ function subfactory_list.add_subfactory(player, name, icon)
 
     local settings = data_util.get("settings", player)
     subfactory.timescale = settings.default_timescale
-    if settings.prefer_matrix_solver then subfactory.matrix_free_items = {} end
+    subfactory.solver_type = settings.default_solver_type
 
     local context = data_util.get("context", player)
     Factory.add(context.factory, subfactory)
diff --git a/modfiles/ui/elements/view_state.lua b/modfiles/ui/elements/view_state.lua
index 8d70c084..2f53a99c 100644
--- a/modfiles/ui/elements/view_state.lua
+++ b/modfiles/ui/elements/view_state.lua
@@ -97,9 +97,9 @@ function view_state.generate_metadata(player, subfactory, formatting_precision,
     }
 end
 
-function view_state.process_item(metadata, item, item_amount, machine_count)
+function view_state.process_item(metadata, item, item_amount, machine_count, is_recipe_line)
     local raw_amount = item_amount or item.amount
-    if raw_amount == nil or (raw_amount < metadata.adjusted_margin_of_error and item.class ~= "Product") then
+    if raw_amount == nil or (raw_amount < metadata.adjusted_margin_of_error and not is_recipe_line) then
         return -1, nil
     end
 
diff --git a/modfiles/ui/event_handler.lua b/modfiles/ui/event_handler.lua
index cd5e39d9..35d7401a 100644
--- a/modfiles/ui/event_handler.lua
+++ b/modfiles/ui/event_handler.lua
@@ -41,6 +41,7 @@ local gui_identifier_map = {
     [defines.events.on_gui_confirmed] = "on_gui_confirmed",
     [defines.events.on_gui_text_changed] = "on_gui_text_changed",
     [defines.events.on_gui_checked_state_changed] = "on_gui_checked_state_changed",
+    [defines.events.on_gui_selection_state_changed] = "on_gui_selection_state_changed",
     [defines.events.on_gui_switch_state_changed] = "on_gui_switch_state_changed",
     [defines.events.on_gui_elem_changed] = "on_gui_elem_changed",
     [defines.events.on_gui_value_changed] = "on_gui_value_changed"
@@ -87,6 +88,10 @@ special_gui_handlers.on_gui_checked_state_changed = (function(event, _, _)
     return true, {state = event.element.state}
 end)
 
+special_gui_handlers.on_gui_selection_state_changed = (function(event, _, _)
+    return true, {selected_index = event.element.selected_index}
+end)
+
 special_gui_handlers.on_gui_switch_state_changed = (function(event, _, _)
     return true, {switch_state = event.element.switch_state}
 end)
diff --git a/modfiles/ui/ui_util.lua b/modfiles/ui/ui_util.lua
index 0244a2de..e278f99b 100755
--- a/modfiles/ui/ui_util.lua
+++ b/modfiles/ui/ui_util.lua
@@ -56,7 +56,7 @@ function ui_util.generate_tutorial_tooltip(player, element_type, has_alt_action,
     local player_table = data_util.get("table", player)
 
     local archive_check = (avoid_archive and player_table.ui_state.flags.archive_open)
-    if player_table.preferences.tutorial_mode and not archive_check then
+    if player_table.preferences.tutorial_mode and not archive_check and element_type then
         local action_tooltip = {"fp.tut_mode_" .. element_type}
 
         local alt_action_name, alt_action_tooltip = player_table.settings.alt_action, ""