(phixonline)-->
with javascript_semantics
function matrix_mul(sequence a, sequence b)
if length(a[1]) != length(b) then
return 0
end if
sequence c = repeat(repeat(0,length(b[1])),length(a))
for i=1 to length(a) do
for j=1 to length(b[1]) do
for k=1 to length(a[1]) do
c[i][j] += a[i][k]*b[k][j]
end for
end for
end for
return c
end function
function pivotize(sequence m)
integer n = length(m)
sequence im = repeat(repeat(0,n),n)
for i=1 to n do
im[i][i] = 1
end for
for i=1 to n do
atom mx = m[i][i]
integer row = i
for j=i to n do
if m[j][i]>mx then
mx = m[j][i]
row = j
end if
end for
if i!=row then
{im[i],im[row]} = {im[row],im[i]}
end if
end for
return im
end function
function lu(sequence a)
integer n = length(a)
sequence l = repeat(repeat(0,n),n),
u = repeat(repeat(0,n),n),
p = pivotize(a),
a2 = matrix_mul(p,a)
for j=1 to n do
l[j][j] = 1.0
for i=1 to j do
atom sum1 = 0.0
for k=1 to i do
sum1 += u[k][j] * l[i][k]
end for
u[i][j] = a2[i][j] - sum1
end for
for i=j+1 to n do
atom sum2 = 0.0
for k=1 to j do
sum2 += u[k][j] * l[i][k]
end for
l[i][j] = (a2[i][j] - sum2) / u[j][j]
end for
end for
return {a, l, u, p}
end function
constant a = {{{1, 3, 5},
{2, 4, 7},
{1, 1, 0}},
{{11, 9,24, 2},
{ 1, 5, 2, 6},
{ 3,17,18, 1},
{ 2, 5, 7, 1}}}
for i=1 to length(a) do
?"== a,l,u,p: =="
pp(lu(a[i]),{pp_Nest,2,pp_Pause,0})
end for