1 files changed, 111 insertions, 0 deletions
diff --git a/lupd.awk b/lupd.awk
new file mode 100644
index 0000000..39bf00c
--- /dev/null
+++ b/lupd.awk
@@ -0,0 +1,111 @@
+#!/usr/bin/awk -f
+
+### lupd.awk
+# LU Decomposition example
+# https://en.wikipedia.org/wiki/LU_decomposition
+# input: square array as delimited text
+# output:
+#
+# note, matrix indexing begins with 1 and uses i=row, j=col
+# lu = L + U - I
+
+# solve A[n,n] * x[n] = b[n]
+function lu_solve(A, b, n) {
+
+    # decomposition of A
+    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"
+        }
+
+    }
+
+    # find solution of Ly = b
+    for (i = 0; i < n; i++) {
+        sum = 0;
+        for (k = 0; k < i; k++)
+            sum += lu[i,k] * y[k];
+        y[i] = b[i] - sum;
+    }
+
+    # find solution of Ux = y
+    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.0 / lu[i,i]) * (y[i] - sum);
+    }
+
+    return x;
+
+}
+
+# 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 (y=1; y<=NF; y++)
+        ($y ~ number) ? header[y] = "col" y : header[y] = $y
+}
+
+# read input data
+NF > 0 {
+    (NF > nf_max) ? nf_max = NF : nf_max = nf_max
+
+    ### columns
+    for (y=1; y<=nf_max; y++) {
+        ### rows
+        if ($y !~ number) {
+            continue
+        }
+        else {
+            read_data[NR,y] = $y
+            col_sum[y] += $y
+            row_sum[NR] += $y
+        }
+    }
+}
+
+END {
+    print length(row_sum), length(col_sum)
+    ### columns
+    for (y=1; y<=nf_max; y++) {
+        if (y in col_sum) {
+            ### rows
+            for (x=1; x<=NR; x++) {
+                if (x in row_sum) {
+                    printf("[" OFMT "," OFMT "]" OFS OFMT OFS OFMT,
+                        x, y, read_data[x,y], col_sum[y])
+                    if (x < nf_max)
+                        printf(OFS)
+                }
+            }
+        }
+        printf(ORS)
+    }
+}
+