Module:LinearAlgebra

local p = {}

function p.MatrixInit(n) local r={} r["mat"]={} for i=1,n do		r["mat"][i]={} end r["m"]=n return r end

--- v_i=1 function p.OnesVector(l) local v={} v["vec"]={} for i=1,l do		v["vec"][i]=1 end v["l"]=l return v end

--- I function p.IdentityMatrix(n) local r=p.MatrixInit(n) r["n"]=n for i = 1,n do		r["mat"][i][i]=1 end return r end

--- R=A*B function p.MatrixMultiplication(A,B) local m1=A["m"] local n1=A["n"] local m2=B["m"] local n2=B["n"] local r = nil if n1==m2 then r=p.MatrixInit(m1) r["m"]=m1 r["n"]=n2 for i = 1,m1 do			for j= 1,n2 do				if A["mat"][i] ~= nil then for k=1,m2 do						if ((A["mat"][i][k] ~= nil) and (B["mat"][k][j] ~= nil)) then if r["mat"][i][j]==nil then r["mat"][i][j]=0 end r["mat"][i][j]= r["mat"][i][j] + A["mat"][i][k] * B["mat"][k][j] end end end end end end return r end

--- r=A*b function p.MatrixVectorMultiplication(A,b) local r = nil if A["n"]==b["l"] then r={} r["vec"]={} r["l"]=A["m"] for i,rows in pairs(A["mat"]) do			for j,entries in pairs(rows) do				r["vec"][i]=0 local z=b["vec"][j] if z ~= nil then r["vec"][i]= r["vec"][i] + entries * z				end end end end return r end

--- R=k*A function p.ScalarMatrixMultiplication(k,A) local r=A for i,rows in pairs(A["mat"]) do		for j,entries in pairs(rows) do			r["mat"][i][j]= k * entries end end return r end

--- R=A+B function p.MatrixSum(A,B) local m1=A["m"] local n1=A["n"] local m2=B["m"] local n2=B["n"] local r = nil if ((n1==n2) and (m1==m2)) then r=A for i=1,m2 do			for j=1,n1 do				local z=B["mat"][i][j] if (z ~= nil) then if r["mat"][i][j] ~= nil then r["mat"][i][j]= r["mat"][i][j] +  z					else r["mat"][i][j]= z					end end end end end return r end

--- R=A-B function p.MatrixDifference(A,B) local m1=A["m"] local n1=A["n"] local m2=B["m"] local n2=B["n"] local r = nil if ((n1==n2) and (m1==m2)) then r=A for i=1,m2 do			for j=1,n1 do				local z=B["mat"][i][j] if (z ~= nil) then if r["mat"][i][j] ~= nil then r["mat"][i][j]= r["mat"][i][j] -  z					else r["mat"][i][j]= -z end end end end end return r end

--- R_ij=A_ij*B_ij function p.MatrixHadamardProduct(A,B) local m1=A["m"] local n1=A["n"] local m2=B["m"] local n2=B["n"] local r = nil if ((n1==n2) and (m1==m2)) then r=p.MatrixInit(m1) r["m"]=m1 r["n"]=n2 for i=1,m2 do			for j=1,n1 do				local u= A["mat"][i][j] local w= B["mat"][i][j] if ((u ~= nil) and (w ~= nil)) then r["mat"][i][j]= u *  w				end end end end return r end

--- r_i=a_i*b_i function p.VectorHadamardProduct(a,b) local l1=a["l"] local l2=b["l"] local r = nil if (l1==l2) then r={} r["vec"]={} r["l"]=l1 for k,v in pairs(a["vec"]) do			local z=b["vec"][k] if z ~= nil then r["vec"][k]= v * b["vec"][k] end end end return r end

--- r=a+b function p.VectorSum(a,b) local l1=a["l"] local l2=b["l"] local r = nil if (l1==l2) then r=a for k,v in pairs(b["vec"]) do			if v ~= nil then if r[k] ~= nil then r[k]= r[k] + v				else r[k] = v				end end end end return r end

--- r=a-b function p.VectorDifference(a,b) local l1=a["l"] local l2=b["l"] local r = nil if (l1==l2) then r=a for k,v in pairs(b["vec"]) do			if v ~= nil then if r[k] ~= nil then r[k]= r[k] - v				else r[k] = -v end end end end return r end

--- adapted from from https://en.wikipedia.org/wiki/LU_decomposition function p.MatrixLUPFactorization(A) local M=#A local N=#A[1] local result = nil if N==M then ---is square local absA local maxA local imax local ptr local Pv={} for i = 1,N do			Pv[i] = i ---Unit permutation matrix, P[N] initialized with N		end

for i = 1,N do			maxA = 0 imax = i

for k = i,N do				absA = math.abs(A[k][i]) if (absA > maxA) then maxA = absA; imax = k;				end end

if (imax ~= i) then ---pivoting P				j = Pv[i] Pv[i] = Pv[imax] Pv[imax] = j

---pivoting rows of A				ptr = A[i] A[i] = A[imax] A[imax] = ptr end

for j = i+1,N do				A[j][i] = A[j][i] / A[i][i]

for k = i+1,N do					A[j][k] = A[j][k] - A[j][i] * A[i][k] end end end

local L=p.IdentityMatrix(N) local U=p.IdentityMatrix(N)

for i=1,N do			for j=1,N do				if j>=i then U[i][j]=A[i][j] else L[i][j]=A[i][j] end end end

local P=p.IdentityMatrix(N) for j = 1,N do			for i = 1,N do				if (Pv[i] == j) then P[i][j] = 1 else P[i][j] = 0 end end end result={} result["L"]=L result["U"]=U result["P"]=P ---decomposition done end return result end

---corrected version from https://algowiki-project.org/en/Forward_substitution --- solve Ly=b function p.ForwardSubstitution(L,B) local lb=#B local M=#L local N=#L[1] local Y = nil if ((lb==M) and (lb==N)) then Y={} Y[1] = B[1] for I = 2,N do			S = B[I] for J = 1,I-1 do				S = S - L[I][J] * Y[J] end Y[I]=S end end return Y end

--- from https://algowiki-project.org/en/Backward_substitution --- solve Ux=y function p.BackwardSubstitution(U,Y) local ly=#Y local M=#U local N=#U[1] local X = nil if ((ly==M) and (ly==N)) then X={} X[N] = Y[N] / U[N][N] for I = N-1,1,-1 do			S = Y[I] for J = N,I+1,-1 do				S = S - U[I][J] * X[J] end X[I] = S / U[I][I] end end return X end

--- solve Ax=b function p.SolveLinearSystem(A,b) local R=p.MatrixLUPFactorization(A) local L=R["L"] local U=R["U"] local P=R["P"] return p.BackwardSubstitution(U,p.ForwardSubstitution(L,p.MatrixVectorMultiplication(1,P,b))) end

function p.LUPInvert(L,U,P) local N=#L local IA=p.MatrixInit(N) local PIv={} local I=p.IdentityMatrix(N) local PI=p.MatrixMultiplication(1,P,I) for j=1,N do		for i=1,N do			PIv[i]=PI[i][j] end local v=p.BackwardSubstitution(U,p.ForwardSubstitution(L,PIv)) for i=1,N do			IA[i][j]=v[i] end end return IA end

--- AI*A=I function p.MatrixInverse(A) local M=#A local N=#A[1] local r = nil if (M==N) then local R=p.MatrixLUPFactorization(A) local L=R["L"] local U=R["U"] local P=R["P"] r = p.LUPInvert(L,U,P) end return r end

function p.MatrixDiagonal(M) local m=#M local n=#M[1] local N=math.min(m,n) local d={} for i=1,N do		d[i]=M[i][i] end return d end

function p.VectorToDiagonalMatrix(v) local n=#v local m=p.IdentityMatrix(n) for i=1,n do		m[i][i]=v[i] end return m end

function p.MatrixRowSwap(M,a,b) local m=#M local n=#M[1] local r = nil if ((a<=m) and (b<=m)) then local R=M for j=1,n do			local t=R[a][j] R[a][j]=R[b][j] R[b][j]=t end r = R	elseif ((a==b) and (a<=m)) then r = M	end return r end

function p.MatrixColumnSwap(M,a,b) local m=#M local n=#M[1] local r = nil if ((a<=n) and (b<=n)) then local R=M for i=1,n do			local t=R[i][a] R[i][a]=R[i][b] R[i][b]=t end r = R	elseif ((a==b) and (a<=n)) then r = M	end return r end

function p.FindInMatrix(M,a) local m=#M local n=#M[1] local r={} local k=0 for i=1,m do		for j=1,n do			if M[i][j]==a then k=k+1 r[k]={i,j} end end end return r end

function p.GetSubmatrix(M,ma,na,mb,nb) local m=#M local n=#M[1] local r = nil if ((na<=nb) and (nb<=n) and (ma<=mb) and (mb<=m)) then local R=p.MatrixInit(mb-ma+1) for i=ma,mb do			for j=na,nb do				R[i-ma+1][j-na+1]=M[i][j] end end r = R	end return r end

function p.MatrixToWikitable(M) local m=M["m"] local n=M["n"] local t='{| class="wikitable"\n' for i=1,m do t = t .. "| "		for j=1,n do t = t .. tostring(M["mat"][i][j]) .. " || " 		end t = t .. "\n|-\n" end t = t .. "\n|}" return t end

return p