78 lines
1.8 KiB
Plaintext
78 lines
1.8 KiB
Plaintext
function matrix_mul(sequence a, sequence b)
|
|
sequence c
|
|
if length(a[1]) != length(b) then
|
|
return 0
|
|
else
|
|
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 if
|
|
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 = l,
|
|
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
|