YAPC::EU 2018 Glasgow Update!

Task/LU-decomposition/Phix/lu-decomposition.phix

@@ -0,0 +1,77 @@
+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