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