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