/*
 * lu.c - LINPACK-style dense solve on obsidian, two ways (C89, gcc 2.95 / xlc).
 * Claude Opus 5.5, 2026-10-11.
 *
 *   lu unblocked N REPS   column-oriented daxpy LU (what LINPACK's dgefa does)
 *   lu blocked N NB       right-looking blocked LU: panel by daxpy, trailing update
 *                         by a 4x4 register-blocked kernel over packed panels
 *
 * Both factor with partial pivoting, solve A x = b with b = A * (1,...,1), check
 * the residual the way LINPACK does, and report MFLOPS for 2n^3/3 + 2n^2 flops.
 * Matrix is column-major, element (i,j) at a[i + j*lda]. Times: gettimeofday,
 * factor+solve only (matrix generation and the copy are outside the timer).
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/time.h>

#define A(i, j) a[(i) + (long)(j) * lda]

static double now(void)
{
    struct timeval tv;
    gettimeofday(&tv, (void *)0);
    return tv.tv_sec + tv.tv_usec / 1e6;
}

static void matgen(double *a, int lda, int n, double *b)
{
    int i, j, init = 1325;
    for (j = 0; j < n; j++)
        for (i = 0; i < n; i++) {
            init = 3125 * init % 65536;
            A(i, j) = (init - 32768.0) / 16384.0;
        }
    for (i = 0; i < n; i++) b[i] = 0;
    for (j = 0; j < n; j++)
        for (i = 0; i < n; i++) b[i] += A(i, j);
}

/* unblocked LU of columns j0..j1-1, rows j0..m-1 (row swaps applied only within
 * these columns; caller swaps the rest) */
static void panel(double *a, int lda, int m, int j0, int j1, int *ipvt)
{
    int i, j, k, p;
    double t, *ck, *cj;
    for (k = j0; k < j1; k++) {
        ck = &A(0, k);
        p = k;
        for (i = k + 1; i < m; i++)
            if (fabs(ck[i]) > fabs(ck[p])) p = i;
        ipvt[k] = p;
        if (p != k)
            for (j = j0; j < j1; j++) { t = A(k, j); A(k, j) = A(p, j); A(p, j) = t; }
        t = -1.0 / ck[k];
        for (i = k + 1; i < m; i++) ck[i] *= t;
        for (j = k + 1; j < j1; j++) {
            cj = &A(0, j);
            t = cj[k];
            if (t != 0)
                for (i = k + 1; i < m; i++) cj[i] += t * ck[i];
        }
    }
}

/* C(m x n) += L(m x kc) * U(kc x n): L packed as 4-row strips (p-major), U as
 * 4-column strips (p-major). m and n multiples of 4 here; edges done by caller. */
#ifndef MC
#define MC 64   /* rows of L kept in the 32 KB L1 while all of U streams past */
#endif
static void kernel(int m, int n, int kc, const double *Lp, const double *Up,
                   double *c, int ldc)
{
    int i, j, p, ic, iend;
    for (ic = 0; ic < m; ic += MC)
    for (j = 0; j < n; j += 4) {
        const double *U0 = Up + (long)j * kc;
        iend = ic + MC < m ? ic + MC : m;
        for (i = ic; i < iend; i += 4) {
            const double *L = Lp + (long)i * kc, *U = U0;
            double c00 = 0, c01 = 0, c02 = 0, c03 = 0, c10 = 0, c11 = 0, c12 = 0, c13 = 0;
            double c20 = 0, c21 = 0, c22 = 0, c23 = 0, c30 = 0, c31 = 0, c32 = 0, c33 = 0;
            double *cc;
            for (p = 0; p < kc; p++) {
                double l0 = L[0], l1 = L[1], l2 = L[2], l3 = L[3];
                double u0 = U[0], u1 = U[1], u2 = U[2], u3 = U[3];
                c00 += l0 * u0; c10 += l1 * u0; c20 += l2 * u0; c30 += l3 * u0;
                c01 += l0 * u1; c11 += l1 * u1; c21 += l2 * u1; c31 += l3 * u1;
                c02 += l0 * u2; c12 += l1 * u2; c22 += l2 * u2; c32 += l3 * u2;
                c03 += l0 * u3; c13 += l1 * u3; c23 += l2 * u3; c33 += l3 * u3;
                L += 4; U += 4;
            }
            cc = c + i + (long)j * ldc;
            cc[0] += c00; cc[1] += c10; cc[2] += c20; cc[3] += c30; cc += ldc;
            cc[0] += c01; cc[1] += c11; cc[2] += c21; cc[3] += c31; cc += ldc;
            cc[0] += c02; cc[1] += c12; cc[2] += c22; cc[3] += c32; cc += ldc;
            cc[0] += c03; cc[1] += c13; cc[2] += c23; cc[3] += c33;
        }
    }
}

static double *Lpack, *Upack;

static void blocked(double *a, int lda, int n, int nb, int *ipvt)
{
    int k, kb, i, j, p, m, nn, m4, n4;
    double t;
    for (k = 0; k < n; k += nb) {
        kb = n - k < nb ? n - k : nb;
        panel(a, lda, n, k, k + kb, ipvt);
        /* apply this panel's row swaps to the columns left and right of it */
        for (p = k; p < k + kb; p++)
            if (ipvt[p] != p) {
                for (j = 0; j < k; j++) { t = A(p, j); A(p, j) = A(ipvt[p], j); A(ipvt[p], j) = t; }
                for (j = k + kb; j < n; j++) { t = A(p, j); A(p, j) = A(ipvt[p], j); A(ipvt[p], j) = t; }
            }
        if (k + kb >= n) break;
        /* U12 = L11^-1 A12 (unit lower; multipliers stored negated, as in dgefa) */
        for (j = k + kb; j < n; j++)
            for (p = k; p < k + kb; p++) {
                t = A(p, j);
                if (t != 0)
                    for (i = p + 1; i < k + kb; i++) A(i, j) += t * A(i, p);
            }
        /* A22 += L21 * U12 (L21 holds negated multipliers) */
        m = n - k - kb; nn = n - k - kb;
        m4 = m & ~3; n4 = nn & ~3;
        for (i = 0; i < m4; i += 4)
            for (p = 0; p < kb; p++) {
                double *d = Lpack + (long)i * kb + p * 4;
                d[0] = A(k + kb + i, k + p); d[1] = A(k + kb + i + 1, k + p);
                d[2] = A(k + kb + i + 2, k + p); d[3] = A(k + kb + i + 3, k + p);
            }
        for (j = 0; j < n4; j += 4)
            for (p = 0; p < kb; p++) {
                double *d = Upack + (long)j * kb + p * 4;
                d[0] = A(k + p, k + kb + j); d[1] = A(k + p, k + kb + j + 1);
                d[2] = A(k + p, k + kb + j + 2); d[3] = A(k + p, k + kb + j + 3);
            }
        kernel(m4, n4, kb, Lpack, Upack, &A(k + kb, k + kb), lda);
        /* ragged edges: rows m4..m-1 of all columns, columns n4..nn-1 of rows < m4 */
        for (j = 0; j < nn; j++)
            for (i = (j < n4 ? m4 : 0); i < m; i++) {
                double s = 0;
                for (p = 0; p < kb; p++) s += A(k + kb + i, k + p) * A(k + p, k + kb + j);
                A(k + kb + i, k + kb + j) += s;
            }
    }
}

/* forward/back substitution like dgesl (job 0), except that row swaps were applied
 * to the whole row (LAPACK style), so all of them are applied to b first */
static void solve(double *a, int lda, int n, int *ipvt, double *b)
{
    int i, k, l;
    double t;
    for (k = 0; k < n - 1; k++) {
        l = ipvt[k];
        if (l != k) { t = b[l]; b[l] = b[k]; b[k] = t; }
    }
    for (k = 0; k < n - 1; k++) {
        t = b[k];
        for (i = k + 1; i < n; i++) b[i] += t * A(i, k);
    }
    for (k = n - 1; k >= 0; k--) {
        b[k] /= A(k, k);
        t = -b[k];
        for (i = 0; i < k; i++) b[i] += t * A(i, k);
    }
}

int main(int argc, char **argv)
{
    int n, lda, nb = 0, reps = 1, r, i, j, *ipvt;
    double *a, *a0, *b, *b0, t, best = 1e30, sum = 0, resid, normx, norma;
    if (argc < 4) { fprintf(stderr, "lu unblocked N REPS | lu blocked N NB\n"); return 2; }
    n = atoi(argv[2]);
    if (!strcmp(argv[1], "blocked")) { nb = atoi(argv[3]); reps = 3; }
    else reps = atoi(argv[3]);
    lda = n + 1;
    a = malloc(sizeof(double) * lda * n); a0 = malloc(sizeof(double) * lda * n);
    b = malloc(sizeof(double) * n); b0 = malloc(sizeof(double) * n);
    ipvt = malloc(sizeof(int) * n);
    Lpack = malloc(sizeof(double) * (n + 4) * (nb + 1));
    Upack = malloc(sizeof(double) * (n + 4) * (nb + 1));
    if (!a || !a0 || !b || !b0 || !ipvt || !Lpack || !Upack) { perror("malloc"); return 1; }
    matgen(a0, lda, n, b0);
    for (r = 0; r < reps; r++) {
        memcpy(a, a0, sizeof(double) * lda * n);
        memcpy(b, b0, sizeof(double) * n);
        t = now();
        if (nb) blocked(a, lda, n, nb, ipvt);
        else panel(a, lda, n, 0, n, ipvt);
        solve(a, lda, n, ipvt, b);
        t = now() - t;
        sum += t;
        if (t < best) best = t;
    }
    /* residual: ||A x - b|| / (n ||A|| ||x|| eps), as LINPACK prints it */
    resid = normx = norma = 0;
    for (i = 0; i < n; i++) {
        double s = -b0[i];
        for (j = 0; j < n; j++) s += a0[i + (long)j * lda] * b[j];
        if (fabs(s) > resid) resid = fabs(s);
        if (fabs(b[i]) > normx) normx = fabs(b[i]);
    }
    for (i = 0; i < n * lda; i++) if (fabs(a0[i]) > norma) norma = fabs(a0[i]);
    t = (2.0 * n * n * n / 3.0 + 2.0 * n * n) / 1e6;
    printf("%s n=%d nb=%d reps=%d: best %.4f s = %.1f MFLOPS, mean %.1f MFLOPS;"
           " norm.resid %.2f, x[0]-1 %.2e\n", argv[1], n, nb, reps, best, t / best,
           t / (sum / reps), resid / (n * norma * normx * 2.220446e-16), b[0] - 1);
    return 0;
}