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

// function initGeMatrix
void
initGeMatrix(size_t m, size_t n,
             double *A,
             ptrdiff_t incRowA, ptrdiff_t incColA)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            A[i*incRowA + j*incColA] = i*n + j +1;
        }
    }
}

// function printGeMatrix
void
printGeMatrix(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");
}

#ifndef COLMAJOR
#define COLMAJOR 1
#endif

int
main()
{
    printf("COLMAJOR = %d\n", COLMAJOR);

    size_t m = 5, n = 10;

    ptrdiff_t incRowA = COLMAJOR ? 1 : n;
    ptrdiff_t incColA = COLMAJOR ? m : 1;

    // allocate memory for A
    double *A = malloc(m*n*sizeof(double));
    if (!A) {
        abort();
    }

    // initialize matrix A
    initGeMatrix(m, n, A, incRowA, incColA);

    // print matrix A
    printf("A =\n");
    printGeMatrix(m, n, A, incRowA, incColA);

    // print how elements are stored in memory
    printf("memory layout of A:\n");
    printGeMatrix(1, m*n, A, 0, 1);

    // print a matrix view of A
    printf("A(0:2, 0:3) =\n");
    printGeMatrix(3, 4, A, incRowA, incColA);

    // print a matrix view of A
    printf("A(2:4, 3:8) =\n");
    printGeMatrix(3, 5, &A[2*incRowA+3*incColA], incRowA, incColA);

    // print a matrix view of A
    printf("A(2:4, 3:8)^T =\n");
    printGeMatrix(5, 3, &A[2*incRowA+3*incColA], incColA, incRowA);

    // release memory
    free(A);
}