#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <math.h>

void
printDGeMatrix(size_t m, size_t n,
               const double *A,
               ptrdiff_t incRowA, ptrdiff_t incColA)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            printf("%9.2lf ", A[i*incRowA + j*incColA]);
        }
        printf("\n");
    }
    printf("\n");
}

void
printIGeMatrix(size_t m, size_t n,
               const size_t *A,
               ptrdiff_t incRowA, ptrdiff_t incColA)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            printf("%9zu ", A[i*incRowA + j*incColA]);
        }
        printf("\n");
    }
    printf("\n");
}

void
dscal(size_t n,
      double alpha,
      double *x, ptrdiff_t incX)
{
    if (alpha==1) {
        return;
    }
    if (alpha!=0) {
        for (size_t i=0; i<n; ++i) {
            x[i*incX] *= alpha;
        }
    } else {
        for (size_t i=0; i<n; ++i) {
            x[i*incX] = 0;
        }
    }
}

void
dger(size_t m, size_t n,
     double alpha,
     const double *x, ptrdiff_t incX,
     const double *y, ptrdiff_t incY,
     double *A, ptrdiff_t incRowA, ptrdiff_t incColA)
{
    if (alpha==0 || m==0 || n==0) {
        return;
    }

    if (incRowA>incColA) {
        dger(n, m, alpha, y, incY, x, incX, A, incColA, incRowA);
        return;
    }

    for (size_t j=0; j<n; ++j) {
        for (size_t i=0; i<m; ++i) {
            A[i*incRowA+j*incColA] += alpha*x[i*incX]*y[j*incY];
        }
    }
}

size_t
idamax(size_t n, const double *x, ptrdiff_t incX)
{
    double max = 0;
    size_t I   = 0;

    for (size_t i=0; i<n; ++i) {
        if (fabs(x[i*incX])>max) {
            I   = i;
            max = fabs(x[i*incX]);
        }
    }
    return I;
}

void
dswap(size_t n, double *x, ptrdiff_t incX, double *y, ptrdiff_t incY)
{
    for (size_t i=0; i<n; ++i) {
        double tmp = x[i*incX];
        x[i*incX] = y[i*incY];
        y[i*incY] = tmp;
    }
}


int
main()
{
    size_t m = 3;
    size_t n = 4;
    size_t min_mn = m<n ? m : n;
    ptrdiff_t incRowA = 1;
    ptrdiff_t incColA = m;
    ptrdiff_t incP = 1;

    double *A = malloc(m*n*sizeof(double));
    size_t *p = malloc(min_mn*sizeof(size_t));
    if (!A || !p) {
        abort();
    }

    A[0*incRowA + 0*incColA] = 2;
    A[0*incRowA + 1*incColA] = 1;
    A[0*incRowA + 2*incColA] = 3;
    A[0*incRowA + 3*incColA] = 1;

    A[1*incRowA + 0*incColA] = 4;
    A[1*incRowA + 1*incColA] = 5;
    A[1*incRowA + 2*incColA] = 6;
    A[1*incRowA + 3*incColA] = 2;

    A[2*incRowA + 0*incColA] = 8;
    A[2*incRowA + 1*incColA] = 8;
    A[2*incRowA + 2*incColA] = 9;
    A[2*incRowA + 3*incColA] = 3;

    printf("A =\n");
    printDGeMatrix(m, n, A, incRowA, incColA);

    // compute LU
    size_t k = m<n ? m : n;

    for (size_t j=0; j<k; ++j) {
        size_t pj = j + idamax(m-j, &A[j*incRowA + j*incColA], incRowA);

        if (j!=pj) {
            dswap(n, &A[j*incRowA], incColA, &A[pj*incRowA], incColA);
        }
        p[j*incP] = pj;

        double alpha = A[j*incRowA + j*incColA];
        if (alpha==0) {
            return j;
        } else {
            alpha = 1/alpha;
        }

        dscal(m-1-j, alpha, &A[(j+1)*incRowA + j*incColA], incRowA);
        dger(m-1-j, n-1-j, -1,
             &A[(j+1)*incRowA +  j   *incColA], incRowA,
             &A[ j   *incRowA + (j+1)*incColA], incColA,
             &A[(j+1)*incRowA + (j+1)*incColA], incRowA, incColA);
    }

    printDGeMatrix(m, n, A, incRowA, incColA);
    printIGeMatrix(1, k, p, 0, incP);

    free(A);
    free(p);
}