1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
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)
}
}
|
