From 309c5d8f7ff8c36235222d079955cd3783bb7ad0 Mon Sep 17 00:00:00 2001
From: wukong <wukong@longaeva>
Date: Mon, 18 Dec 2023 12:37:31 -0800
Subject: Added a partially working function performing LU decomposition on a
 square matrix, ludcmp.awk. This is providing correct answers on _some_ test
 data, further testing debugging required here.

---
 ludcmp.awk | 105 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 105 insertions(+)
 create mode 100644 ludcmp.awk

(limited to 'ludcmp.awk')

diff --git a/ludcmp.awk b/ludcmp.awk
new file mode 100644
index 0000000..6c7807a
--- /dev/null
+++ b/ludcmp.awk
@@ -0,0 +1,105 @@
+#!/usr/bin/awk -f
+
+### ludcmp.awk
+# LU decomposition on square matrix, based on example code from,
+# https://en.wikipedia.org/wiki/LU_decomposition
+# INPUT:
+#   A (N-by-N matrix)
+#   right_part (N-element vector)
+# OUTPUT:
+#   x (N-element vector)
+# note, matrix indexing begins with 1 and uses i=row, j=col
+# lu = L + U - I
+
+# decomposition of A
+function lu_decomp(A, n) {
+
+    for (i=1; i<=n; i++) {
+
+        for (j=i; j<=n; j++) {
+            sum = 0.0
+            for (k=1; k<=i; k++)
+                sum += lu[i,k] * lu[k,j]
+            lu[i,j] = A[i,j] - sum
+        }
+
+        for (j=i+1; j<=n; j++) {
+            sum = 0.0
+            for (k=1; k<=i; k++)
+                sum += lu[j,k] * lu[k,i]
+
+            # check for div by zero
+            if (lu[i,i] != 0.0)
+                lu[j,i] = (A[j,i] - sum) / lu[i,i]
+            else
+                lu[j,i] = "nan"
+        }
+
+        # debug
+        for (j=1; j<=n; j++) {
+            printf(lu[i,j])
+            printf( j==n ? ORS : OFS )
+        }
+    }
+}
+
+# regex to identify strings that look like numbers
+BEGIN {
+    OFS = FS
+    sign = "[+-]?"
+    decimal = "[0-9]+[.]?[0-9]*"
+    fraction = "[.][0-9]*"
+    exponent = "([Ee]" sign "[0-9]+)?"
+    number = "^" sign "(" decimal "|" fraction ")" exponent "$"
+}
+
+# column headers
+NR == 1 {
+    for (n=1; n<=NF; n++)
+        ($n ~ number) ? head[n] = "col" n : head[n] = $n
+}
+
+# read input data
+NF {
+
+    if (NF > max_nf)
+        (NF > nf_max) ? nf_max = NF : nf_max = nf_max
+
+    for (n=1; n<=NF; n++) {
+        if ($n ~ number) {
+            count[n] += 1
+            (count[n] == 1 || max_nf > size) ? size = max_nf : size = size
+            (count[n] == 1 || count[n] > size) ? size = count[n] : size = size
+            matrix[count[n],n] = $n
+            print(matrix[count[n],n])
+        }
+    }
+
+}
+
+END {
+    printf(ORS)
+    print(NR, max_nf, size)
+    printf(ORS)
+    lu_decomp(matrix, size)
+
+}
+
+## find solution of Ly = b, for y
+#for (i=0; i<n;, i++) {
+#    sum = 0
+#    for (k=0; k<i; k++)
+#        sum += lu[i,k] * y[k]
+#    y[i] = right_part[i] - sum
+#}
+#
+## find solution of Ux = y, for x
+#for (i=n-1; i>=0; i--) {
+#    sum = 0;
+#    for (k=i+1; k<n; k++)
+#        sum += lu[i,k] * x[k];
+#    x[i] = (1 / lu[i,i]) * (y[i] - sum);
+#}
+#
+#return x
+
-- 
cgit v1.2.3