50 lines
1.2 KiB
Text
50 lines
1.2 KiB
Text
F pprint(m)
|
||
L(row) m
|
||
print(row)
|
||
|
||
F matrix_mul(a, b)
|
||
V result = [[0.0] * a.len] * a.len
|
||
L(j) 0 .< a.len
|
||
L(i) 0 .< a.len
|
||
V r = 0.0
|
||
L(k) 0 .< a.len
|
||
r += a[i][k] * b[k][j]
|
||
result[i][j] = r
|
||
R result
|
||
|
||
F pivotize(m)
|
||
‘Creates the pivoting matrix for m.’
|
||
V n = m.len
|
||
V ID = (0 .< n).map(j -> (0 .< @n).map(i -> Float(i == @j)))
|
||
L(j) 0 .< n
|
||
V row = max(j .< n, key' i -> abs(@m[i][@j]))
|
||
I j != row
|
||
swap(&ID[j], &ID[row])
|
||
R ID
|
||
|
||
F lu(A)
|
||
‘Decomposes a nxn matrix A by PA=lU and returns l, U and P.’
|
||
V n = A.len
|
||
V l = [[0.0] * n] * n
|
||
V U = [[0.0] * n] * n
|
||
V P = pivotize(A)
|
||
V A2 = matrix_mul(P, A)
|
||
L(j) 0 .< n
|
||
l[j][j] = 1.0
|
||
L(i) 0 .. j
|
||
V s1 = sum((0 .< i).map(k -> @U[k][@j] * @l[@i][k]))
|
||
U[i][j] = A2[i][j] - s1
|
||
L(i) j .< n
|
||
V s2 = sum((0 .< j).map(k -> @U[k][@j] * @l[@i][k]))
|
||
l[i][j] = (A2[i][j] - s2) / U[j][j]
|
||
R (l, U, P)
|
||
|
||
V a = [[1, 3, 5], [2, 4, 7], [1, 1, 0]]
|
||
L(part) lu(a)
|
||
pprint(part)
|
||
print()
|
||
print()
|
||
V b = [[11, 9, 24, 2], [1, 5, 2, 6], [3, 17, 18, 1], [2, 5, 7, 1]]
|
||
L(part) lu(b)
|
||
pprint(part)
|
||
print()
|