local M={}local g={}local b=function(M)local b=#M local h=#M[1]local t=setmetatable({},g)for g=1,b,1 do local b=h*(g-1)for h=1,h,1 do if M[g][h]~=0 then t[b+h]=M[g][h]end end end rawset(t,"\100\105\109\101\110\115\105\111\110\115",{b,h})return t end local h=function(M,b)b=b or 0 local h=#M local t=#M[1]local V=setmetatable({},g)for g=1,h,1 do local h=t*(g-1)for t=1,t,1 do if math.abs(M[g][t])>=b then V[h+t]=M[g][t]end end end rawset(V,"\100\105\109\101\110\115\105\111\110\115",{h;t})return V end local t=function(M)local b=setmetatable({},g)rawset(b,"\100\105\109\101\110\115\105\111\110\115",{M.dimensions[1],M.dimensions[2]})for M,g in pairs(M)do if type(M)=="\110\117\109\098\101\114"then b[M]=g end end return b end local V=function(M)local b=setmetatable({},g)rawset(b,"\100\105\109\101\110\115\105\111\110\115",{M,M})for g=1,M,1 do b[M*(g-1)+g]=1 end return b end local k=function(M,b)local h=setmetatable({},g)rawset(h,"\100\105\109\101\110\115\105\111\110\115",{M,b})return h end local B=function(M,b,h,t,V)V=V or 0 local k=setmetatable({},g)rawset(k,"\100\105\109\101\110\115\105\111\110\115",{M;b})for M=1,M,1 do local g=b*(M-1)for M=1,b,1 do local b=(t-h)*math.random()+h if math.abs(b)>=V then k[g+M]=b end end end return k end local Q=function(M,g,b)local h=M.dimensions[1]local t=M.dimensions[2]for h=1,h,1 do local V=t*(h-1)M[V+g],M[V+b]=M[V+b],M[V+g]end return M end local d=function(M,g)local b=V(M)for M=1,M,1 do local h=g[M]if h~=nil then b=Q(b,M,h)end end return b end local n=function(M,g,b)local h=M.dimensions[1]local t=M.dimensions[2]for h=1,t,1 do local V=t*(g-1)local k=t*(b-1)M[V+h],M[k+h]=M[k+h],M[V+h]end return M end local L=function(M,g)local b=V(M)for M=1,M,1 do local h=g[M]if h~=nil then b=n(b,M,h)end end return b end local f=function(M,g,b,h,t)t=t or 0 local V=M.dimensions[2]for k=1,V,1 do local B=M[V*(g-1)+k]local Q=M[V*(b-1)+k]if B~=nil and Q~=nil then local g=Q+h*B if math.abs(g)<t then M[V*(b-1)+k]=nil else M[V*(b-1)+k]=g end elseif B~=nil then local g=h*B if math.abs(g)<t then M[V*(b-1)+k]=nil else M[V*(b-1)+k]=g end end end return M end local G=function(M,g,b,h,t)t=t or 0 local V=M.dimensions[2]for k=1,V,1 do local B=M[V*(b-1)+k]local Q=g[k]if B~=nil and Q~=nil then local g=B+h*Q if math.abs(g)<t then M[V*(b-1)+k]=nil else M[V*(b-1)+k]=g end elseif Q~=nil then local g=h*Q if math.abs(g)<t then M[V*(b-1)+k]=nil else M[V*(b-1)+k]=g end end end return M end local e=function(M)local b=M.dimensions[1]local h=M.dimensions[2]local t=setmetatable({},g)rawset(t,"\100\105\109\101\110\115\105\111\110\115",{h,b})for M,g in pairs(M)do if type(M)=="\110\117\109\098\101\114"then local V=M%h if V==0 then V=h end local k=(M-V)/h+1 t[b*(V-1)+k]=g end end return t end local c=function(M)local g=M.dimensions[1]local b=M.dimensions[2]if g~=b then error("\067\097\110\110\111\116 \099\111\109\112\117\116\101 \076\085 \111\102 \115\112\097\114\115\101 \114\101\099\116\097\110\103\117\108\097\114 \109\097\116\114\105\120\046")end local h={}local t=V(g)for V=1,b-1,1 do local k=V local B=M[b*(V-1)+V]or 0 for g=V,g,1 do local h=M[b*(g-1)+V]or 0 h=math.abs(h)if h~=nil and h>B then B=h k=g end end if B==0 then error("\083\112\097\114\115\101 \109\097\116\114\105\120 \105\115 \110\111\116 \105\110\118\101\114\116\105\098\108\101")end if k~=V then n(M,V,k)for M=1,V-1,1 do t[b*(V-1)+M],t[b*(k-1)+M]=t[b*(k-1)+M],t[b*(V-1)+M]end h[V]=k end B=M[b*(V-1)+V]for g=V+1,g,1 do local h=M[b*(g-1)+V]if h~=nil then local k=h/B t[b*(g-1)+V]=k M[b*(g-1)+V]=nil for h=V+1,b,1 do local t=M[b*(g-1)+h]local B=M[b*(V-1)+h]if t~=nil and B~=nil then M[b*(g-1)+h]=t-B*k elseif B~=nil then M[b*(g-1)+h]=-B*k end end end end end local k=L(g,h)return{t;M,k}end local Z=function(M)local g=M.dimensions[1]local b=M.dimensions[2]if g~=b then error("\067\097\110\110\111\116 \099\111\109\112\117\116\101 \105\110\118\101\114\115\101 \111\102 \115\112\097\114\115\101 \114\101\099\116\097\110\103\117\108\097\114 \109\097\116\114\105\120\046")end local h=V(g)for t=1,b-1,1 do local V=t local k=M[b*(t-1)+t]or 0 for g=t,g,1 do local h=M[b*(g-1)+t]or 0 h=math.abs(h)if h~=nil and h>k then k=h V=g end end if k==0 then error("\083\112\097\114\115\101 \109\097\116\114\105\120 \105\115 \110\111\116 \105\110\118\101\114\116\105\098\108\101")end if V~=t then n(M,t,V)n(h,t,V)end k=M[b*(t-1)+t]for g=t+1,g,1 do local V=M[b*(g-1)+t]if V~=nil then local B=V/k M[b*(g-1)+t]=nil f(h,t,g,-B)for h=t+1,b,1 do local V=M[b*(g-1)+h]local k=M[b*(t-1)+h]if V~=nil and k~=nil then M[b*(g-1)+h]=V-k*B elseif k~=nil then M[b*(g-1)+h]=-k*B end end end end end for g=g,1,-1 do local t=b*(g-1)local V=M[t+g]for t=1,g-1,1 do local k=M[b*((g-t)-1)+g]if k~=nil then f(h,g,g-t,-k/V)end end for M=1,b,1 do h[t+M]=h[t+M]/V end end return h end local W=function(M,g)g=g or 0 for b,h in pairs(M)do if type(b)=="\110\117\109\098\101\114"then if h<g then M[b]=nil end end end return M end local Y={}local X=function(M)local g=setmetatable({},Y)rawset(g,"\107\101\121\115",{})rawset(g,"\108\101\110\103\116\104",0)for b=1,#M,1 do if M[b]~=0 then g[b]=M[b]g.keys[g.keys.length+1]=b end end return g end local C=function(M,g)g=g or 0 local b=setmetatable({},Y)rawset(b,"\107\101\121\115",{})rawset(b,"\108\101\110\103\116\104",#M)for h=1,#M,1 do if math.abs(M[h])>=g then b[h]=M[h]b.keys[#b.keys+1]=h end end return b end local O=function(M)local g=setmetatable({},Y)rawset(g,"\107\101\121\115",{})rawset(g,"\108\101\110\103\116\104",M.length)for M in M.keys do g.keys[#g.keys+1]=M end return g end local T=function(M,g)local b=setmetatable({},Y)rawset(b,"\107\101\121\115",{g})rawset(b,"\108\101\110\103\116\104",M)b[g]=1 return b end local r=function(M)local g=setmetatable({},Y)rawset(g,"\107\101\121\115",{})rawset(g,"\108\101\110\103\116\104",M)return g end local A=function(M,g,b,h)h=h or 0 local t=setmetatable({},Y)rawset(t,"\107\101\121\115",{})rawset(t,"\108\101\110\103\116\104",M)for M=1,M,1 do local V=(b-g)*math.random()+g if V>=h then t[M]=V t.keys[#t.keys+1]=M end end return t end local u=function(M)local g=setmetatable({},Y)rawset(g,"\107\101\121\115",{})rawset(g,"\108\101\110\103\116\104",M)for M=1,M,1 do g[M]=1 g.keys[#g.keys+1]=M end return g end local P=function(M,g)local b=1 if#M<#g then M,g=g,M end local h={}local t=1 for V=1,#M,1 do local k=M[V]local B=g[b]while B and k>B do h[t]=B t=t+1 b=b+1 B=g[b]end if k==B then h[t]=k t=t+1 b=b+1 else h[t]=k t=t+1 end end return h end local s=function(M,g,b)M[g],M[b]=M[b],M[g]return M end local a={}local I=function(M)local g=setmetatable({},a)rawset(g,"\100\105\109\101\110\115\105\111\110\115",{#M;#M[1]})rawset(g,"\107\101\121\115",{})for b=1,#M,1 do local h=X(M[b])if#h.keys~=0 then g[b]=h g.keys[#g.keys+1]=b end end return g end local H=function(M)local g=setmetatable({},a)rawset(g,"\100\105\109\101\110\115\105\111\110\115",{M;M})rawset(g,"\107\101\121\115",{})for b=1,M,1 do g[b]=T(M,b)g.keys[#g.keys+1]=b end return g end local m=function(M,g)for g,b in ipairs(g)do M[g],M[b]=M[b],M[g]end local b={}for M,g in ipairs(M)do b[#b+1]=M end M.keys=b return M end local v=function(M,g)local b=H(M)return m(b,g)end local F=function(M) end local K=function(M,g)local b=setmetatable({},Y)local h={}rawset(b,"\108\101\110\103\116\104",M.dimensions[1])for t,V in ipairs(M.keys)do b[V]=M[V][g]h[#h+1]=V end rawset(b,"\107\101\121\115",h)return b end local E=function(M)local g={}for b,h in ipairs(M.keys)do g=P(g,M[h].keys)end return g end local p=function(M)local b=setmetatable({},g)local h=M.dimensions[1]local t=M.dimensions[2]rawset(b,"\100\105\109\101\110\115\105\111\110\115",{h,t})for g,h in ipairs(M.keys)do local V=M[h]local k=t*(h-1)for M,g in ipairs(V.keys)do b[k+g]=V[g]end end return b end local z=function(M)local g=setmetatable({},a)local b=M.dimensions[1]local h=M.dimensions[2]rawset(g,"\100\105\109\101\110\115\105\111\110\115",{b,h})local t={}local V=0 for M,b in pairs(M)do if type(M)=="\110\117\109\098\101\114"then local k=M%h if k==0 then k=h end local B=math.floor((M-k)/h+1.1)if B>V then local M=r(h)M[k]=b M.keys[#M.keys+1]=k g[B]=M t[#t+1]=B V=B else local M=g[V]M[k]=b M.keys[#M.keys+1]=k t[#t+1]=V end end end rawset(g,"\107\101\121\115",t)return g end local y={}local R=function(M)local g=setmetatable(M,y)rawset(g,"\100\105\109\101\110\115\105\111\110\115",{#M;#M[1]})return g end local N=function(M,g)local b={}for g=1,g,1 do if g==M then b[g]=1 else b[g]=0 end end return b end local j=function(M)local g={}for M=1,M,1 do g[M]=0 end return g end local D=function(M)local g={}for b=1,M,1 do g[b]=N(b,M)end return R(g)end local l=function(M,g)local b={}for M=1,M,1 do b[M]=j(g)end return R(b)end Y.__add=function(M,g)if M.length~=g.length then error("\067\097\110\110\111\116 \097\100\100 \115\112\097\114\115\101 \118\101\099\116\111\114\115 \111\102 \117\110\101\113\117\097\108 \108\101\110\103\116\104\046",2)end local b=P(M.keys,g.keys)local h=setmetatable({},Y)rawset(h,"\108\101\110\103\116\104",M.length)local t={}for b,V in ipairs(b)do local k=M[V]local B=g[V]if k and B then local M=k+B if M==0 then h[V]=nil else h[V]=M t[#t+1]=V end elseif B then h[V]=B t[#t+1]=V else h[V]=k t[#t+1]=V end end h.keys=t return h end Y.__sub=function(M,g)if M.length~=g.length then error("\067\097\110\110\111\116 \097\100\100 \115\112\097\114\115\101 \118\101\099\116\111\114\115 \111\102 \117\110\101\113\117\097\108 \108\101\110\103\116\104\046",2)end local b=P(M.keys,g.keys)local h=setmetatable({},Y)rawset(h,"\108\101\110\103\116\104",M.length)local t={}for b,V in ipairs(b)do local k=M[V]local B=g[V]if k and B then local M=k-B if M==0 then h[V]=nil else h[V]=M t[#t+1]=V end elseif B then h[V]=-B t[#t+1]=V else h[V]=k t[#t+1]=V end end h.keys=t return h end Y.__unm=function(M)local g=O(M)for M,b in ipairs(g.keys)do g[b]=-g[b]end return g end Y.__mul=function(M,g)if M.length~=g.length then error("\067\097\110\110\111\116 \116\097\107\101 \105\110\110\101\114 \112\114\111\100\117\099\116 \111\102 \115\112\097\114\115\101 \118\101\099\116\111\114\115 \111\102 \117\110\101\113\117\097\108 \108\101\110\103\116\104\046",2)end local b=P(M.keys,g.keys)local h=0 for b,t in ipairs(b)do local V=M[t]local k=g[t]if V and k then h=h+V*k end end return h end Y.__tostring=function(M)local g="\123"local b=M.length for b=1,b-1,1 do local h=M[b]if h then g=g..(h.."\044")else g=g.."\048\044"end end local h=M[b]if h then g=g..h else g=g..0 end g=g.."\125"return g end a.__add=function(M,g)if M.dimensions[1]~=g.dimensions[1]and M.dimensions[2]~=g.dimensions[2]then error("\067\097\110\110\111\116 \097\100\100 \115\112\097\114\115\101 \109\097\116\114\105\099\101\115 \111\102 \117\110\101\113\117\097\108 \100\105\109\101\110\115\105\111\110\115\046",2)end local b=P(M.keys,g.keys)local h=setmetatable({},a)rawset(h,"\100\105\109\101\110\115\105\111\110\115",{M.dimensions[1];M.dimensions[2]})local t={}for b,V in ipairs(b)do local k=M[V]local B=g[V]if k and B then local M=k+B if M==0 then h[V]=nil else h[V]=M t[#t+1]=V end elseif B then h[V]=B t[#t+1]=V else h[V]=k t[#t+1]=V end end h.keys=t return h end a.__sub=function(M,g)if M.dimensions[1]~=g.dimensions[1]and M.dimensions[2]~=g.dimensions[2]then error("\067\097\110\110\111\116 \097\100\100 \115\112\097\114\115\101 \109\097\116\114\105\099\101\115 \111\102 \117\110\101\113\117\097\108 \100\105\109\101\110\115\105\111\110\115\046",2)end local b=P(M.keys,g.keys)local h=setmetatable({},a)rawset(h,"\100\105\109\101\110\115\105\111\110\115",{M.dimensions[1];M.dimensions[2]})local t={}for b,V in ipairs(b)do local k=M[V]local B=g[V]if k and B then local M=k-B if M==0 then h[V]=nil else h[V]=M t[#t+1]=V end elseif B then h[V]=-B t[#t+1]=V else h[V]=k t[#t+1]=V end end h.keys=t return h end a.__mul=function(M,g)local b=p(M)local h=p(g)local t=z(b*h)return t end a.__tostring=function(M)local g="\123"local b=M.dimensions[1]local h=M.dimensions[2]for b=1,b-1,1 do local t=M[b]if t then g=g..(tostring(t).."\044")else g=g..(tostring(r(h)).."\044")end end local t=M[b]if t then g=g..tostring(t)else g=g..tostring(r(h))end g=g.."\125"return g end g.__add=function(M,g)if M.dimensions[1]~=g.dimensions[1]or M.dimensions[2]~=g.dimensions[2]then error("\065\116\116\101\109\112\116\105\110\103 \116\111 \097\100\100 \115\112\097\114\115\101 \109\097\116\114\105\099\101\115 \111\102 \100\105\102\102\101\114\101\110\116 \115\105\122\101\115\046")end local b=z(M)local h=z(g)return p(b+h)end g.__sub=function(M,g)if M.dimensions[1]~=g.dimensions[1]or M.dimensions[2]~=g.dimensions[2]then error("\065\116\116\101\109\112\116\105\110\103 \116\111 \097\100\100 \115\112\097\114\115\101 \109\097\116\114\105\099\101\115 \111\102 \100\105\102\102\101\114\101\110\116 \115\105\122\101\115\046")end local b=z(M)local h=z(g)return p(b-h)end g.__mul=function(M,b)if M.dimensions[2]~=b.dimensions[1]then error("\065\116\116\101\109\112\116\105\110\103 \116\111 \109\117\108\116\105\112\108\121 \105\110\099\111\109\112\097\116\105\098\108\101 \115\112\097\114\115\101 \109\097\116\114\105\099\101\115\046")end local h=setmetatable({},g)rawset(h,"\100\105\109\101\110\115\105\111\110\115",{M.dimensions[1];b.dimensions[2]})for g,t in pairs(M)do if type(g)=="\110\117\109\098\101\114"then local V=g%M.dimensions[2]if V==0 then V=M.dimensions[2]end local k=(g-V)/M.dimensions[2]local B=b.dimensions[2]*(V-1)local Q=b.dimensions[2]*k for M=1,b.dimensions[2],1 do local g=h[Q+M]local V=b[B+M]if type(g)=="\110\117\109\098\101\114"and V~=nil then g=g+t*b[B+M]h[Q+M]=g elseif V~=nil then h[Q+M]=t*b[B+M]end end end end return h end g.__tostring=function(M)local g="\123"local b=M.dimensions[1]local h=M.dimensions[2]for b=1,b-1,1 do g=g.."\123"for t=1,h-1,1 do local V=M[h*(b-1)+t]or 0 g=g..(tostring(V).."\044")end local t=M[h*(b-1)+h]or 0 g=g..(tostring(t).."\125\044")end g=g.."\123"for t=1,h-1,1 do local V=M[h*(b-1)+t]or 0 g=g..(tostring(V).."\044")end local t=M[h*(b-1)+h]or 0 g=g..(tostring(t).."\125")g=g.."\125"return g end y.__add=function(M,g)if M.dimensions[1]~=g.dimensions[1]or M.dimensions[2]~=g.dimensions[2]then error("\065\116\116\101\109\112\116\105\110\103 \116\111 \097\100\100 \109\097\116\114\105\099\101\115 \111\102 \100\105\102\102\101\114\101\110\116 \115\105\122\101\115\046")end local b={}for M,h in ipairs(M)do b[M]={}for h,t in ipairs(h)do b[M]=t+g[M][h]end end return R(b)end y.__sub=function(M,g)if M.dimensions[1]~=g.dimensions[1]or M.dimensions[2]~=g.dimensions[2]then error("\065\116\116\101\109\112\116\105\110\103 \116\111 \097\100\100 \109\097\116\114\105\099\101\115 \111\102 \100\105\102\102\101\114\101\110\116 \115\105\122\101\115\046")end local b={}for M,h in ipairs(M)do b[M]={}for h,t in ipairs(h)do b[M]=t-g[M][h]end end return R(b)end M.sparseVector={}M.sparseVector.new=function(M)return X(M)end M.sparseVector.newNumeric=function(M,g)return C(M,g)end M.sparseVector.ones=function(M)return u(M)end M.sparseVector.zero=function(M)return r(M)end M.sparseVector.random=function(M,g,b,h)return A(M,g,b,h)end M.sparseMatrix={}M.sparseMatrix.new=function(M)return I(M)end M.sparseMatrix.identity=function(M)return H(M)end M.sparseMatrix.flatten=function(M)return p(M)end M.liSparseMatrix={}M.liSparseMatrix.new=function(M)return b(M)end M.liSparseMatrix.newNumeric=function(M,g)return h(M,g)end M.liSparseMatrix.identity=function(M)return V(M)end M.liSparseMatrix.zero=function(M,g)return k(M,g)end M.liSparseMatrix.random=function(M,g,b,h,t)return B(M,g,b,h,t)end M.liSparseMatrix.lu=function(M)local g=t(M)return c(g)end M.liSparseMatrix.inverse=function(M)local g=t(M)return Z(g)end M.liSparseMatrix.transpose=function(M)local g=t(M)return e(g)end M.liSparseMatrix.sparsify=function(M,g)local b=t(M)return W(b,g)end M.liSparseMatrix.unflatten=function(M)return z(M)end M.matrix={}M.matrix.new=function(M)return R(M)end M.matrix.identity=function(M)return D(M)end M.matrix.zero=function(M,g)return l(M,g)end return M
local matrixAlgebra = {}

local _linearlyIndexedSparseMatrix = {}

local _linearlyIndexedSparseMatrixFromTableOfTables = function (tableOfTables)
    local numberOfRows = #tableOfTables
    local numberOfColumns = #tableOfTables[1]

    local sparseForm = setmetatable({}, _linearlyIndexedSparseMatrix)

    for i = 1, numberOfRows do
        local rowConstant = numberOfColumns * (i - 1)
        for j = 1, numberOfColumns do
            if tableOfTables[i][j] ~= 0 then
                sparseForm[rowConstant + j] = tableOfTables[i][j]
            end
        end
    end

    rawset(sparseForm, "dimensions", {numberOfRows, numberOfColumns})

    return sparseForm
end

local _linearlyIndexedSparseMatrixFromNumericTableOfTables = function (tableOfTables, tol)
    tol = tol or 0

    local numberOfRows = #tableOfTables
    local numberOfColumns = #tableOfTables[1]

    local sparseForm = setmetatable({}, _linearlyIndexedSparseMatrix)

    for i = 1, numberOfRows do
        local rowConstant = numberOfColumns * (i - 1)
        for j = 1, numberOfColumns do
            if math.abs(tableOfTables[i][j]) >= tol then
                sparseForm[rowConstant + j] = tableOfTables[i][j]
            end
        end
    end

    rawset(sparseForm, "dimensions", {numberOfRows, numberOfColumns})

    return sparseForm
end

local _linearlyIndexedSparseCopy = function (matrix)
    local copy = setmetatable({}, _linearlyIndexedSparseMatrix)

    rawset(copy, "dimensions", {matrix.dimensions[1], matrix.dimensions[2]})

    for k, v in pairs(matrix) do
        if type(k) == "number" then
            copy[k] = v
        end
    end

    return copy
end

local _linearlyIndexedSparseIdentity = function (n)
    local result = setmetatable({}, _linearlyIndexedSparseMatrix)

    rawset(result, "dimensions", {n, n})

    for i = 1, n do
        result[n * (i - 1) + i] = 1
    end

    return result
end

local _linearlyIndexedSparseZero = function (n, m)
    local result = setmetatable({}, _linearlyIndexedSparseMatrix)

    rawset(result, "dimensions", {n,m})

    return result
end

local _linearlyIndexedSparseRandom = function (n, m, a, b, tol)
    tol = tol or 0

    local result = setmetatable({}, _linearlyIndexedSparseMatrix)

    rawset(result, "dimensions", {n, m})

    for i = 1, n do
        local rowConstant = m * (i - 1)
        for j = 1, m do
            local val = (b - a) * (math.random()) + a
            if math.abs(val) >= tol then
                result[rowConstant + j] = val
            end
        end
    end

    return result
end

local _linearlyIndexedSparseColumnSwap = function (matrix, n, m)
    local numberOfRows = matrix.dimensions[1]
    local numberOfColumns = matrix.dimensions[2]

    for i = 1, numberOfRows do
        local rowConstant = numberOfColumns * (i - 1)
        matrix[rowConstant + n],  matrix[rowConstant + m] =  matrix[rowConstant + m],  matrix[rowConstant + n]
    end

    return matrix
end

local _linearlyIndexedSparseRightPermutationMatrix = function (n, permutations)
    local result = _linearlyIndexedSparseIdentity(n)

    for i = 1, n do
        local destination = permutations[i]
        if destination ~= nil then
            result = _linearlyIndexedSparseColumnSwap(result, i, destination)
        end
    end

    return result
end

local _linearlyIndexedSparseRowSwap = function (matrix, n, m)
    local numberOfRows = matrix.dimensions[1]
    local numberOfColumns = matrix.dimensions[2]

    for i = 1, numberOfColumns do
        local rowConstant = numberOfColumns * (n - 1)
        local rowConstant2 = numberOfColumns * (m - 1)
        matrix[rowConstant + i],  matrix[rowConstant2 + i] =  matrix[rowConstant2 + i],  matrix[rowConstant + i]
    end

    return matrix
end

local _linearlyIndexedSparsePermutationMatrix = function (n, permutations)
    local result = _linearlyIndexedSparseIdentity(n)

    for i = 1, n do
        local destination = permutations[i]
        if destination ~= nil then
            result = _linearlyIndexedSparseRowSwap(result, i, destination)
        end
    end

    return result
end

local _linearlyIndexedSparseShear = function (matrix, n, m, c, tol)
    tol = tol or 0

    local numberOfColumns = matrix.dimensions[2]

    for i = 1, numberOfColumns do
        local val1 = matrix[numberOfColumns * (n - 1) + i]
        local val2 = matrix[numberOfColumns * (m - 1) + i]
        if val1 ~= nil and val2 ~= nil then
            local val3 = val2 + c * val1
            if math.abs(val3) < tol then
                matrix[numberOfColumns * (m - 1) + i] = nil
            else
                matrix[numberOfColumns * (m - 1) + i] = val3
            end
        elseif val1 ~= nil then
            local val3 = c * val1
            if math.abs(val3) < tol then
                matrix[numberOfColumns * (m - 1) + i] = nil
            else
                matrix[numberOfColumns * (m - 1) + i] = val3
            end
        end
    end

    return matrix
end

local _linearlyIndexedSparseRowAdd = function (matrix, row, n, c, tol)
    tol = tol or 0

    local numberOfColumns = matrix.dimensions[2]

    for i = 1, numberOfColumns do
        local val1 = matrix[numberOfColumns * (n - 1) + i]
        local val2 = row[i]
        if val1 ~= nil and val2 ~= nil then
            local val3 = val1 + c * val2
            if math.abs(val3) < tol then
                matrix[numberOfColumns * (n - 1) + i] = nil
            else
                matrix[numberOfColumns * (n - 1) + i] = val3
            end
        elseif val2 ~= nil then
            local val3 = c * val2
            if math.abs(val3) < tol then
                matrix[numberOfColumns * (n - 1) + i] = nil
            else
                matrix[numberOfColumns * (n - 1) + i] = val3
            end
        end
    end

    return matrix
end

local _linearlyIndexedSparseTranspose = function (matrix)
    local numberOfRows = matrix.dimensions[1]
    local numberOfColumns = matrix.dimensions[2]

    local result = setmetatable({}, _linearlyIndexedSparseMatrix)

    rawset(result, "dimensions", {numberOfColumns, numberOfRows})

    for k, v in pairs(matrix) do
        if type(k) == "number" then
            local columnNumber = k % numberOfColumns
            if columnNumber == 0 then
                columnNumber = numberOfColumns
            end
            local rowNumber = (k - columnNumber) / numberOfColumns + 1
            result[numberOfRows * (columnNumber - 1) + rowNumber] = v
        end
    end

    return result
end

local _linearlyIndexedSparseLU = function (matrix)
    local numberOfRows = matrix.dimensions[1]
    local numberOfColumns = matrix.dimensions[2]

    if numberOfRows ~= numberOfColumns then
        error("Cannot compute LU of sparse rectangular matrix.")
    end

    local permuations = {}

    local l = _linearlyIndexedSparseIdentity(numberOfRows)

    for i = 1, numberOfColumns - 1 do
        local maxRow = i
        local max = matrix[numberOfColumns * (i - 1) + i] or 0

        for j = i, numberOfRows do
            local maxCandidate = matrix[numberOfColumns * (j - 1) + i] or 0
            maxCandidate = math.abs(maxCandidate)
            if maxCandidate ~= nil and maxCandidate > max then
                max = maxCandidate
                maxRow = j
            end
        end

        if max == 0 then
            error("Sparse matrix is not invertible")
        end

        if maxRow ~= i then
            _linearlyIndexedSparseRowSwap(matrix, i, maxRow)
            for k = 1, i - 1 do
                l[numberOfColumns * (i - 1) + k], l[numberOfColumns * (maxRow - 1) + k] = l[numberOfColumns * (maxRow - 1) + k], l[numberOfColumns * (i - 1) + k]
            end
            permuations[i] = maxRow
        end

        max = matrix[numberOfColumns * (i - 1) + i]

        for j = i + 1, numberOfRows do
            local val = matrix[numberOfColumns * (j - 1) + i]
            if val ~= nil then
                local valOverMax = val / max
                l[numberOfColumns * (j - 1) + i] = valOverMax
                matrix[numberOfColumns * (j - 1) + i] = nil
                for k = i + 1, numberOfColumns do
                    local val1 = matrix[numberOfColumns * (j - 1) + k]
                    local val2 = matrix[numberOfColumns * (i - 1) + k]
                    if val1 ~= nil and val2 ~= nil then
                        matrix[numberOfColumns * (j - 1) + k] = val1 - val2 * valOverMax
                    elseif val2 ~= nil then
                        matrix[numberOfColumns * (j - 1) + k] = -val2 * valOverMax
                    end
                end
            end
        end
    end

    local permuationMatrix = _linearlyIndexedSparsePermutationMatrix(numberOfRows, permuations)

    return {l, matrix, permuationMatrix}
end

local _linearlyIndexedSparseInverse = function (matrix)
    local numberOfRows = matrix.dimensions[1]
    local numberOfColumns = matrix.dimensions[2]

    if numberOfRows ~= numberOfColumns then
        error("Cannot compute inverse of sparse rectangular matrix.")
    end

    local result = _linearlyIndexedSparseIdentity(numberOfRows)

    for i = 1, numberOfColumns - 1 do
        local maxRow = i
        local max = matrix[numberOfColumns * (i - 1) + i] or 0

        for j = i, numberOfRows do
            local maxCandidate = matrix[numberOfColumns * (j - 1) + i] or 0
            maxCandidate = math.abs(maxCandidate)
            if maxCandidate ~= nil and maxCandidate > max then
                max = maxCandidate
                maxRow = j
            end
        end

        if max == 0 then
            error("Sparse matrix is not invertible")
        end

        if maxRow ~= i then
            _linearlyIndexedSparseRowSwap(matrix, i, maxRow)
            _linearlyIndexedSparseRowSwap(result, i, maxRow)
        end

        max = matrix[numberOfColumns * (i - 1) + i]

        for j = i + 1, numberOfRows do
            local val = matrix[numberOfColumns * (j - 1) + i]
            if val ~= nil then
                local valOverMax = val / max
                matrix[numberOfColumns * (j - 1) + i] = nil
                _linearlyIndexedSparseShear(result, i, j, -valOverMax)
                for k = i + 1, numberOfColumns do
                    local val1 = matrix[numberOfColumns * (j - 1) + k]
                    local val2 = matrix[numberOfColumns * (i - 1) + k]
                    if val1 ~= nil and val2 ~= nil then
                        matrix[numberOfColumns * (j - 1) + k] = val1 - val2 * valOverMax
                    elseif val2 ~= nil then
                        matrix[numberOfColumns * (j - 1) + k] = -val2 * valOverMax
                    end
                end
            end
        end
    end

    for i = numberOfRows, 1, -1 do
        local rowConstant = numberOfColumns * (i - 1)
        local val = matrix[rowConstant + i]
        for j = 1, i - 1 do
            local val1 = matrix[numberOfColumns * (i - j - 1) + i]
            if val1 ~= nil then
                _linearlyIndexedSparseShear(result, i, i - j, -val1 / val)
            end
        end
        for j = 1, numberOfColumns do
            result[rowConstant + j] = result[rowConstant + j] / val
        end
    end

    return result
end

local _linearlyIndexedSparseMatrixSparsify = function (matrix, tol)
    tol = tol or 0

    for k, v in pairs(matrix) do
        if type(k) == "number" then
            if v < tol then
                matrix[k] = nil
            end
        end
    end

    return matrix
end

local _sparseVector = {}

local _sparseVectorFromArray = function (array)
    local result = setmetatable({}, _sparseVector)

    rawset(result, "keys", {})

    rawset(result, "length", 0)

    for i = 1, #array do
        if array[i] ~= 0 then
            result[i] = array[i]
            result.keys[result.keys.length + 1] = i
        end
    end

    return result
end

local _sparseVectorFromNumericArray = function (array, tol)
    tol = tol or 0

    local result = setmetatable({}, _sparseVector)

    rawset(result, "keys", {})

    rawset(result, "length", #array)

    for i = 1, #array do
        if math.abs(array[i]) >= tol then
            result[i] = array[i]
            result.keys[#result.keys + 1] = i
        end
    end

    return result
end

local _sparseVectorCopy = function (vector)
    local copy = setmetatable({}, _sparseVector)

    rawset(copy, "keys", {})

    rawset(copy, "length", vector.length)

    for k in vector.keys do
        copy.keys[#copy.keys + 1] = k
    end

    return copy
end

local _sparseStandardBasisVector = function (n, m)
    local result = setmetatable({}, _sparseVector)

    rawset(result, "keys", {m})

    rawset(result, "length", n)

    result[m] = 1

    return result
end

local _sparseZeroVector = function (n)
    local result = setmetatable({}, _sparseVector)

    rawset(result, "keys", {})

    rawset(result, "length", n)

    return result
end

local _sparseRandomVector = function (n, a, b, tol)
    tol = tol or 0

    local result = setmetatable({}, _sparseVector)

    rawset(result, "keys", {})

    rawset(result, "length", n)

    for i = 1, n do
        local val = (b - a) * (math.random()) + a
        if val >= tol then
            result[i] = val
            result.keys[#result.keys + 1] = i
        end
    end

    return result
end

local _sparseOnesVector = function (n)
    local result = setmetatable({}, _sparseVector)

    rawset(result, "keys", {})

    rawset(result, "length", n)

    for i = 1, n do
        result[i] = 1
        result.keys[#result.keys + 1] = i
    end

    return result
end

local _keyUnion = function (left, right)
    local j = 1

    if #left < #right then
        left, right = right, left
    end

    local keys = {}

    local index = 1

    for i = 1, #left do
        local val1 = left[i]
        local val2 = right[j]

        while val2 and val1 > val2 do
            keys[index] = val2
            index = index + 1
            j = j + 1
            val2 = right[j]
        end
        if val1 == val2 then
            keys[index] = val1
            index = index + 1
            j = j + 1
        else
            keys[index] = val1
            index = index + 1
        end
    end

    return keys
end

local _sparseVectorSwap = function (vector, n, m)
    vector[n], vector[m] = vector[m], vector[n]

    return vector
end

local _sparseMatrix = {}

local _sparseMatrixFromArrayOfArrays = function (ArrayOfArrays)
    local result = setmetatable({}, _sparseMatrix)

    rawset(result, "dimensions", {#ArrayOfArrays, #ArrayOfArrays[1]})

    rawset(result, "keys", {})

    for i = 1, #ArrayOfArrays do
        local row = _sparseVectorFromArray(ArrayOfArrays[i])
        if #row.keys ~= 0 then
            result[i] = row
            result.keys[#result.keys+1] = i
        end
    end

    return result
end

local _sparseMatrixIdentity = function (n)
    local result = setmetatable({}, _sparseMatrix)

    rawset(result, "dimensions", {n, n})

    rawset(result, "keys", {})

    for i = 1, n do
        result[i] = _sparseStandardBasisVector(n, i)
        result.keys[#result.keys+1] = i
    end

    return result
end

local _sparseMatrixRowPermute = function (matrix, permuations)
    for n, m in ipairs(permuations) do
        matrix[n], matrix[m] = matrix[m], matrix[n]
    end

    local keys = {}

    for k, v in ipairs(matrix) do
        keys[#keys+1] = k
    end

    matrix.keys = keys

    return matrix
end

local _sparseLeftPermutationMatrix = function (n, permuations)
    local result = _sparseMatrixIdentity(n)

    return _sparseMatrixRowPermute(result, permuations)
end

local _sparseMatrixLU = function (matrix)
    
end

local _sparseMatrixColumn = function (matrix, n)
    local result = setmetatable({}, _sparseVector)

    local keys = {}

    rawset(result, "length", matrix.dimensions[1])

    for k, i in ipairs(matrix.keys) do
        result[i] = matrix[i][n]
        keys[#keys+1] = i
    end

    rawset(result, "keys", keys)

    return result
end

local _sparseMatrixColumnKeys = function (matrix)
    local keys = {}

    for k, i in ipairs(matrix.keys) do
        keys = _keyUnion(keys, matrix[i].keys)
    end

    return keys
end

local _sparseMatrixFlatten = function (matrix)
    local result = setmetatable({}, _linearlyIndexedSparseMatrix)

    local numberOfRows = matrix.dimensions[1]
    local numberOfColumns = matrix.dimensions[2]

    rawset(result, "dimensions", {numberOfRows, numberOfColumns})

    for k, i in ipairs(matrix.keys) do
        local row = matrix[i]
        local rowConstant = numberOfColumns * (i - 1)
        for kk, ii in ipairs(row.keys) do
            result[rowConstant + ii] = row[ii]
        end
    end

    return result
end

local _sparseMatrixUnFlatten = function (matrix)
    local result = setmetatable({}, _sparseMatrix)

    local numberOfRows = matrix.dimensions[1]
    local numberOfColumns = matrix.dimensions[2]

    rawset(result, "dimensions", {numberOfRows, numberOfColumns})

    local keys = {}

    local rowNumber = 0

    for k, v in pairs(matrix) do
        if type(k) == "number" then
            local columnNumber = k % numberOfColumns
            if columnNumber == 0 then
                columnNumber = numberOfColumns
            end
            local tempRowNumber = math.floor((k - columnNumber) / numberOfColumns + 1.1)
            if tempRowNumber > rowNumber then
                local row = _sparseZeroVector(numberOfColumns)
                row[columnNumber] = v
                row.keys[#row.keys+1] = columnNumber
                result[tempRowNumber] = row
                keys[#keys+1] = tempRowNumber
                rowNumber = tempRowNumber
            else
                local row = result[rowNumber]
                row[columnNumber] = v
                row.keys[#row.keys+1] = columnNumber
                keys[#keys+1] = rowNumber
            end
        end
    end

    rawset(result, "keys", keys)

    return result
end

local _matrix = {}

local _matrixFromTableOfTables = function (tableOfTables)
    local result = setmetatable(tableOfTables, _matrix)

    rawset(result, "dimensions", {#tableOfTables,#tableOfTables[1]})

    return result
end

local _nthStandardBasisVector = function (n, m)
    local e = {}

    for i = 1, m do
        if i == n then
            e[i] = 1
        else
            e[i] = 0
        end
    end

    return e
end

local _zeroVector = function (n)
    local z = {}

    for i = 1, n do
        z[i] = 0
    end

    return z
end

local _matrixIdentity = function (n)
    local result = {}

    for i = 1, n do
        result[i] = _nthStandardBasisVector(i, n)
    end

    return _matrixFromTableOfTables(result)
end

local _matrixZero = function (n, m)
    local result = {}

    for i = 1, n do
        result[i] = _zeroVector(m)
    end

    return _matrixFromTableOfTables(result)
end

_sparseVector.__add = function (left, right)
    if left.length ~= right.length then
        error("Cannot add sparse vectors of unequal length.", 2)
    end

    local keys = _keyUnion(left.keys, right.keys)

    local result = setmetatable({}, _sparseVector)

    rawset(result, "length", left.length)

    local newKeys = {}

    for k, i in ipairs(keys) do
        local val1 = left[i]
        local val2 = right[i]
        if val1 and val2 then
            local val3 = val1 + val2
            if val3 == 0 then
                result[i] = nil
            else
                result[i] = val3
                newKeys[#newKeys + 1] = i
            end
        elseif val2 then
            result[i] = val2
            newKeys[#newKeys + 1] = i
        else
            result[i] = val1
            newKeys[#newKeys + 1] = i
        end
    end

    result.keys = newKeys

    return result
end

_sparseVector.__sub = function (left, right)
    if left.length ~= right.length then
        error("Cannot add sparse vectors of unequal length.", 2)
    end

    local keys = _keyUnion(left.keys, right.keys)

    local result = setmetatable({}, _sparseVector)

    rawset(result, "length", left.length)

    local newKeys = {}

    for k, i in ipairs(keys) do
        local val1 = left[i]
        local val2 = right[i]
        if val1 and val2 then
            local val3 = val1 - val2
            if val3 == 0 then
                result[i] = nil
            else
                result[i] = val3
                newKeys[#newKeys + 1] = i
            end
        elseif val2 then
            result[i] = -val2
            newKeys[#newKeys + 1] = i
        else
            result[i] = val1
            newKeys[#newKeys + 1] = i
        end
    end

    result.keys = newKeys

    return result
end

_sparseVector.__unm = function (vector)
    local copy = _sparseVectorCopy(vector)

    for k, i in ipairs(copy.keys) do
        copy[i] = -copy[i]
    end

    return copy
end

_sparseVector.__mul = function (left, right)
    if left.length ~= right.length then
        error("Cannot take inner product of sparse vectors of unequal length.", 2)
    end

    local keys = _keyUnion(left.keys, right.keys)

    local result = 0

    for k, i in ipairs(keys) do
        local val1 = left[i]
        local val2 = right[i]
        if val1 and val2 then
            result = result + val1 * val2
        end
    end

    return result
end

_sparseVector.__tostring = function (vector)
    local result = "{"

    local length = vector.length

    for i = 1, length - 1 do
        local val = vector[i]
        if val then
            result = result .. val .. ","
        else
            result = result .. 0 .. ","
        end
    end

    local val = vector[length]
    if val then
        result = result .. val
    else
        result = result .. 0
    end

    result = result .. "}"

    return result
end

_sparseMatrix.__add = function (left, right)
    if left.dimensions[1] ~= right.dimensions[1] and left.dimensions[2] ~= right.dimensions[2] then
        error("Cannot add sparse matrices of unequal dimensions.", 2)
    end

    local keys = _keyUnion(left.keys, right.keys)

    local result = setmetatable({}, _sparseMatrix)

    rawset(result, "dimensions", {left.dimensions[1], left.dimensions[2]})

    local newKeys = {}

    for k, i in ipairs(keys) do
        local val1 = left[i]
        local val2 = right[i]
        if val1 and val2 then
            local val3 = val1 + val2
            if val3 == 0 then
                result[i] = nil
            else
                result[i] = val3
                newKeys[#newKeys + 1] = i
            end
        elseif val2 then
            result[i] = val2
            newKeys[#newKeys + 1] = i
        else
            result[i] = val1
            newKeys[#newKeys + 1] = i
        end
    end

    result.keys = newKeys

    return result
end

_sparseMatrix.__sub = function (left, right)
    if left.dimensions[1] ~= right.dimensions[1] and left.dimensions[2] ~= right.dimensions[2] then
        error("Cannot add sparse matrices of unequal dimensions.", 2)
    end

    local keys = _keyUnion(left.keys, right.keys)

    local result = setmetatable({}, _sparseMatrix)

    rawset(result, "dimensions", {left.dimensions[1], left.dimensions[2]})

    local newKeys = {}

    for k, i in ipairs(keys) do
        local val1 = left[i]
        local val2 = right[i]
        if val1 and val2 then
            local val3 = val1 - val2
            if val3 == 0 then
                result[i] = nil
            else
                result[i] = val3
                newKeys[#newKeys + 1] = i
            end
        elseif val2 then
            result[i] = -val2
            newKeys[#newKeys + 1] = i
        else
            result[i] = val1
            newKeys[#newKeys + 1] = i
        end
    end

    result.keys = newKeys

    return result
end

_sparseMatrix.__mul = function (left, right)
    local liLeft = _sparseMatrixFlatten(left)
    local liRight = _sparseMatrixFlatten(right)

    local result = _sparseMatrixUnFlatten(liLeft * liRight)

    return result
end

_sparseMatrix.__tostring = function (matrix)
    local result = "{"

    local length = matrix.dimensions[1]
    local width = matrix.dimensions[2]

    for i = 1, length - 1 do
        local val = matrix[i]
        if val then
            result = result .. tostring(val) .. ","
        else
            result = result .. tostring(_sparseZeroVector(width)) .. ","
        end
    end

    local val = matrix[length]
    if val then
        result = result .. tostring(val)
    else
        result = result .. tostring(_sparseZeroVector(width))
    end

    result = result .. "}"

    return result
end

_linearlyIndexedSparseMatrix.__add = function (left, right)
    if left.dimensions[1] ~= right.dimensions[1] or left.dimensions[2] ~= right.dimensions[2] then
        error("Attempting to add sparse matrices of different sizes.")
    end

    local inflatedLeft = _sparseMatrixUnFlatten(left)
    local inflatedRight = _sparseMatrixUnFlatten(right)

    return _sparseMatrixFlatten(inflatedLeft + inflatedRight)
end

_linearlyIndexedSparseMatrix.__sub = function (left, right)
    if left.dimensions[1] ~= right.dimensions[1] or left.dimensions[2] ~= right.dimensions[2] then
        error("Attempting to add sparse matrices of different sizes.")
    end

    local inflatedLeft = _sparseMatrixUnFlatten(left)
    local inflatedRight = _sparseMatrixUnFlatten(right)

    return _sparseMatrixFlatten(inflatedLeft - inflatedRight)
end

_linearlyIndexedSparseMatrix.__mul = function (left, right)
    if left.dimensions[2] ~= right.dimensions[1] then
        error("Attempting to multiply incompatible sparse matrices.")
    end

    local result = setmetatable({}, _linearlyIndexedSparseMatrix)

    rawset(result, "dimensions", {left.dimensions[1], right.dimensions[2]})

    for k, v in pairs(left) do
        if type(k) == "number" then
            local columnNumber = k % left.dimensions[2]
            if columnNumber == 0 then
                columnNumber = left.dimensions[2]
            end
            local rowNumber = (k - columnNumber) / left.dimensions[2]
            local rowConstant = right.dimensions[2] * (columnNumber - 1)
            local rowConstant2 = right.dimensions[2] * rowNumber
            for i = 1, right.dimensions[2] do
                local val1 = result[rowConstant2 + i]
                local val2 = right[rowConstant + i]
                if type(val1) == "number" and val2 ~= nil then
                    val1 = val1 + v * right[rowConstant + i]
                    result[rowConstant2 + i] = val1
                elseif val2 ~= nil then
                    result[rowConstant2 + i] = v * right[rowConstant + i]
                end
            end
        end
    end

    return result
end

_linearlyIndexedSparseMatrix.__tostring = function (matrix)
    local result = "{"

    local numberOfRows = matrix.dimensions[1]
    local numberOfColumns = matrix.dimensions[2]

    for i = 1, numberOfRows - 1 do
        result = result .. "{"
        for j = 1, numberOfColumns - 1 do
            local val = matrix[numberOfColumns * (i - 1) + j] or 0
            result = result .. tostring(val) .. ","
        end
        local val = matrix[numberOfColumns * (i - 1) + numberOfColumns] or 0
        result = result .. tostring(val) .. "},"
    end

    result = result .. "{"
    for j = 1, numberOfColumns - 1 do
        local val = matrix[numberOfColumns * (numberOfRows - 1) + j] or 0
        result = result .. tostring(val) .. ","
    end
    local val = matrix[numberOfColumns * (numberOfRows - 1) + numberOfColumns] or 0
    result = result .. tostring(val) .. "}"

    result = result .. "}"

    return result
end

_matrix.__add = function (left, right)
    if left.dimensions[1] ~= right.dimensions[1] or left.dimensions[2] ~= right.dimensions[2] then
        error("Attempting to add matrices of different sizes.")
    end

    local result = {}

    for k, v in ipairs(left) do
        result[k] = {}
        for kk, vv in ipairs(v) do
            result[k] = vv + right[k][kk]
        end
    end

    return _matrixFromTableOfTables(result)
end

_matrix.__sub = function (left, right)
    if left.dimensions[1] ~= right.dimensions[1] or left.dimensions[2] ~= right.dimensions[2] then
        error("Attempting to add matrices of different sizes.")
    end

    local result = {}

    for k, v in ipairs(left) do
        result[k] = {}
        for kk, vv in ipairs(v) do
            result[k] = vv - right[k][kk]
        end
    end

    return _matrixFromTableOfTables(result)
end

matrixAlgebra.sparseVector = {}

matrixAlgebra.sparseVector.new = function (array)
    return _sparseVectorFromArray(array)
end

matrixAlgebra.sparseVector.newNumeric = function (array, tol)
    return _sparseVectorFromNumericArray(array, tol)
end

matrixAlgebra.sparseVector.ones = function (n)
    return _sparseOnesVector(n)
end

matrixAlgebra.sparseVector.zero = function (n)
    return _sparseZeroVector(n)
end

matrixAlgebra.sparseVector.random = function (n, a, b, tol)
    return _sparseRandomVector(n, a, b, tol)
end

matrixAlgebra.sparseMatrix = {}

matrixAlgebra.sparseMatrix.new = function (ArrayOfArrays)
    return _sparseMatrixFromArrayOfArrays(ArrayOfArrays)
end

matrixAlgebra.sparseMatrix.identity = function (n)
    return _sparseMatrixIdentity(n)
end

matrixAlgebra.sparseMatrix.flatten = function (matrix)
    return _sparseMatrixFlatten(matrix)
end

matrixAlgebra.liSparseMatrix = {}

matrixAlgebra.liSparseMatrix.new = function (tableOfTables)
    return _linearlyIndexedSparseMatrixFromTableOfTables(tableOfTables)
end

matrixAlgebra.liSparseMatrix.newNumeric = function (tableOfTables, tol)
    return _linearlyIndexedSparseMatrixFromNumericTableOfTables(tableOfTables, tol)
end

matrixAlgebra.liSparseMatrix.identity = function (n)
    return _linearlyIndexedSparseIdentity(n)
end

matrixAlgebra.liSparseMatrix.zero = function (n, m)
    return _linearlyIndexedSparseZero(n, m)
end

matrixAlgebra.liSparseMatrix.random = function (n, m, a, b, tol)
    return _linearlyIndexedSparseRandom(n, m, a, b, tol)
end

matrixAlgebra.liSparseMatrix.lu = function (input)
    local matrix = _linearlyIndexedSparseCopy(input)
    return _linearlyIndexedSparseLU(matrix)
end

matrixAlgebra.liSparseMatrix.inverse = function (input)
    local matrix = _linearlyIndexedSparseCopy(input)
    return _linearlyIndexedSparseInverse(matrix)
end

matrixAlgebra.liSparseMatrix.transpose = function (input)
    local matrix = _linearlyIndexedSparseCopy(input)
    return _linearlyIndexedSparseTranspose(matrix)
end

matrixAlgebra.liSparseMatrix.sparsify = function (input, tol)
    local matrix = _linearlyIndexedSparseCopy(input)
    return _linearlyIndexedSparseMatrixSparsify(matrix, tol)
end

matrixAlgebra.liSparseMatrix.unflatten = function (matrix)
    return _sparseMatrixUnFlatten(matrix)
end

matrixAlgebra.matrix = {}

matrixAlgebra.matrix.new = function (tableOfTables)
    return _matrixFromTableOfTables(tableOfTables)
end

matrixAlgebra.matrix.identity = function (n)
    return _matrixIdentity(n)
end

matrixAlgebra.matrix.zero = function (n, m)
    return _matrixZero(n, m)
end

return matrixAlgebra
