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()