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
|
#!/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
|
