#include <stdio.h>
#include <stdlib.h>
#include <stddef.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
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];
        }
    }
}

int
main()
{
    size_t m = 3;
    size_t n = 4;
    ptrdiff_t incRowA = 1;
    ptrdiff_t incColA = m;

    double *A = malloc(m*n*sizeof(double));
    if (!A) {
        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) {
        double alpha = 1/A[j*incRowA + j*incColA];

        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);
    }

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

    free(A);
}