typedef struct nml_mat_s {
  unsigned int num_rows;
  unsigned int num_cols;
  double **data;
  int is_square;
} nml_mat;

nml_mat *m = ...
m->data[i][j] = ...

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {


  nml_mat* m, *mx;

  printf("\nCreating an empty matrix with 2x3\n");
  m = nml_mat_new(2,3);
  nml_mat_print(m);
  nml_mat_free(m);

  printf("\nCreating a square matrix 5x5 \n");
  m = nml_mat_sqr(5);
  nml_mat_print(m);
  nml_mat_free(m);

  printf("\nCreating an ID 7x7 Matrix and copying it into another matrix:\n");
  m = nml_mat_eye(7);
  mx = nml_mat_cp(m);
  nml_mat_print(m);
  nml_mat_print(mx);
  nml_mat_free(m);
  nml_mat_free(mx);

  return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) { 
    double array[6] = { 
        1.0, 0.2, 3.0, 4.0, 5.0, 3.1 
    };
    nml_mat* my;
    
    // 3 rows, 2 columns
    // read exactly 6 numbers from array[6]
    my = nml_mat_from(3, 2, 6, array);
    nml_mat_print(my);
    nml_mat_free(my);

    // 4 rows, 2 columns
    // read exactly 3 numbers from array[6]
    my = nml_mat_from(4, 2, 3, array);
    nml_mat_print(my);
    nml_mat_free(my);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    const char *f = "examples/data/matrix1.data";
    nml_mat *from_file = nml_mat_fromfile(f);
    nml_mat_print(from_file);
    nml_mat_free(from_file);

    // Or if the file is already opened

    FILE *m_file = fopen("examples/data/matrix2.data", "r");
    nml_mat *from_file2 = nml_mat_fromfilef(m_file);
    nml_mat_print(from_file2);
    nml_mat_free(from_file2);
    fclose(m_file);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *from_file2 = nml_mat_fromfilef(stdin);
    nml_mat_print(from_file2);
    nml_mat_free(from_file2);
    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
  srand(time(NULL)); // Should be called once per program
  nml_mat *m = nml_mat_rnd(5, 5, -10.0, 10.0);
  nml_mat_print(m);
  nml_mat_free(m);
  return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {

    srand(time(NULL));
    nml_mat *m1 = nml_mat_rnd(2, 3, 1.0, 10.0);
    nml_mat *m2 = nml_mat_rnd(2, 3, 1.0, 10.0);

    if (nml_mat_eq(m1, m2, 0.001)) {
        printf("Wow, what were the oddss..\n");
    } else {
        printf("It's ok, nobody is that lucky!\n");
    }
    if (nml_mat_eqdim(m1, m2)) {
        printf("At least we know they both have the same number of rows and columns.\n");
    }

    nml_mat_free(m1);
    nml_mat_free(m2);
    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
  printf("\nExtract all matrix columns from a Matrix as matrices\n");
  srand(time(NULL));
  nml_mat *m = nml_mat_rnd(5, 5, -10.0, 10.0);
  nml_mat *col;
  nml_mat_print(m);
  int i = 0;
  for(i = 0; i < m->num_cols; i++) {
    col = nml_mat_col_get(m, i);
    nml_mat_print(col);
    nml_mat_free(col);
  }
  nml_mat_free(m);
  return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    // Creates a matrix of zeros of size = 5
    nml_mat *pi_mat = nml_mat_sqr(5);

    // Sets all elements to PI
    nml_mat_all_set(pi_mat, M_PI);

    nml_mat_print(pi_mat);
    nml_mat_free(pi_mat);
    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    // Creates a matrix of zeros of size = 5
    nml_mat *pi_mat = nml_mat_sqr(5);

    // Sets the first diagonal to PI
    nml_mat_diag_set(pi_mat, M_PI);

    nml_mat_print(pi_mat);
    nml_mat_free(pi_mat);
    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *m = nml_mat_eye(5);

    // Multiply all elements of m with 2.0
    // and return a new matrix

    nml_mat *new_m = nml_mat_smult(m, 2.0);

    if (!(nml_mat_eq(m, new_m, 0.0))) {
        printf("It's normal to see this message.\n");
    }

    // Multiply all elements of m with 2.0
    // m is modified, no new matrix is created
    nml_mat_smult_r(m, 2.0);

    if (nml_mat_eq(m, new_m, 0.0)) {
        printf("It's even more normal to see this message.\n");
    }

    nml_mat_free(m);
    nml_mat_free(new_m);
    
    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *a = nml_mat_new(4,5);
    nml_mat_all_set(a, 1.0);
    int i = 0;
    for(; i < a->num_rows; ++i) {
        // Changes are doing on matrix a
        // row[i] is multiplied with (double) i
        nml_mat_row_mult_r(a, i, (double)i);
    }
    nml_mat_print(a);

    // Create a new matrix b by multiplying row[1] 
    // in matrix a with 5.0.
    // Matrix a remains unchanged
    nml_mat *b = nml_mat_row_mult(a, 1, 5.0);
    nml_mat_print(b);
    nml_mat_free(a);
    nml_mat_free(b);
    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {

    nml_mat *m = nml_mat_rnd(5, 4, 1.0, 2.0);
    nml_mat_print(m);

    // Add row[1] elements to row[2] elements

    nml_mat_row_addrow_r(m, 2, 1, 1.0);

    // Add row[1] to row[0] with a multiplier of 2.0

    nml_mat_row_addrow_r(m, 0, 1, 2.0);

    nml_mat_print(m);
    nml_mat_free(m);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {

    nml_mat *m = nml_mat_sqr_rnd(4, 1.0, 2.0);
    nml_mat_print(m);

    // Remove column[1] from m
    // m remains the same
    // less_columns is another matrix 
    nml_mat *less_columns = nml_mat_col_rem(m, 1);
    nml_mat_print(less_columns);

    // Remove row[0] from less_columns
    // less_columns remains the same
    // less_rows is another matrix
    nml_mat *less_rows = nml_mat_row_rem(less_columns, 0);
    nml_mat_print(less_rows);

    nml_mat_free(m);
    nml_mat_free(less_columns);
    nml_mat_free(less_rows);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
  srand(time(NULL));

  nml_mat *m = nml_mat_sqr_rnd(8, 0.0, 10.0);

  printf("m=");
  nml_mat_print(m);

  printf("m= (...after swapping col1=%d with col2=%d):\n", 1, 2);
  nml_mat_col_swap_r(m, 1, 2);
  nml_mat_print(m);

  printf("newm= (...after swapping col1=%d with col2=%d and creating a new matrix):\n", 0, 1);
  nml_mat *newm = nml_mat_col_swap(m, 0, 1);
  nml_mat_print(newm);

  printf("m= (...after swapping row1=%d with row2=%d)\n", 0, 2);
  nml_mat_row_swap_r(m, 0, 2);
  nml_mat_print(m);

  nml_mat_free(m);
  nml_mat_free(newm);

  return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
  nml_mat *I = nml_mat_eye(3);
  nml_mat *Ix2 = nml_mat_smult(I, 2.0);
  nml_mat *rndm = nml_mat_rnd(3, 4, 1.0, 5.0);

  nml_mat **ms = malloc(sizeof(*ms) * 2);
  ms[0] = I;
  ms[1] = Ix2;
  
  nml_mat *concats1 = nml_mat_cath(2, ms);

  ms[0] = concats1;
  ms[1] = rndm;

  nml_mat *concats2 = nml_mat_cath(2, ms);

  printf("\nConcatenate horizontally\n");
  printf("I=\n");
  nml_mat_print(I);
  printf("Ix2=\n");
  nml_mat_print(Ix2);
  printf("rndm=\n");
  nml_mat_print(rndm);
  printf("concats1=\n");
  nml_mat_print(concats1);
  printf("concats2=\n");
  nml_mat_print(concats2);

  free(ms);
  nml_mat_free(I);
  nml_mat_free(Ix2);
  nml_mat_free(concats1);
  nml_mat_free(concats2);
  nml_mat_free(rndm);

  // -------------------------------------
  // Vertical concatenation
  // -------------------------------------

  nml_mat *A = nml_mat_rnd(3, 4, 1.0, 4.0);
  nml_mat *B = nml_mat_rnd(5, 4, 10.0, 20.0);
  nml_mat *C = nml_mat_eye(4);

  nml_mat **ABarr = malloc(sizeof(*ABarr) * 2);
  ABarr[0] = A;
  ABarr[1] = B;
  nml_mat *ABCat = nml_mat_catv(2, ABarr);

  printf("\nA=\n");
  nml_mat_print(A);
  printf("\nB=\n");
  nml_mat_print(B);
  printf("\nC=\n");
  nml_mat_print(C);
  printf("\nA concat B =\n");
  nml_mat_print(ABCat);

  free(ABarr);
  nml_mat_free(A);
  nml_mat_free(B);
  nml_mat_free(C);

  return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *m1 = nml_mat_sqr_rnd(4, 0.0, 10.0);
    nml_mat *m2 = nml_mat_sqr_rnd(4, 0.0, 10.0);
    printf("m1=\n");
    nml_mat_print(m1);
    
    printf("m2=\n");
    nml_mat_print(m2);


    // Add the matrices to, result is kept in m3
    // m1 and m2 remain unchanged
    nml_mat *m3 = nml_mat_add(m1, m2);
    printf("m3=\n");
    nml_mat_print(m3);

    // Add the matrices, the result is kept in m1
    // m1 is modified, m2 remains unchanged
    nml_mat_add_r(m1, m2);
    printf("m1=\n");
    nml_mat_print(m1);
    
    nml_mat_free(m1);
    nml_mat_free(m2);
    nml_mat_free(m3);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *m1 = nml_mat_sqr_rnd(4, 0.0, 10.0);
    nml_mat *m2 = nml_mat_sqr_rnd(4, 0.0, 10.0);

    printf("m1=\n");
    nml_mat_print(m1);
    
    printf("m2=\n");
    nml_mat_print(m2);

    // Multiply matrices
    nml_mat *m3 = nml_mat_dot(m1, m2);
    printf("m3=\n");
    nml_mat_print(m3);

    nml_mat_free(m1);
    nml_mat_free(m2);
    nml_mat_free(m3);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *m1 = nml_mat_rnd(1, 5, 1.0, 10.0);
    nml_mat_print(m1);

    nml_mat *m2 = nml_mat_transp(m1);
    nml_mat_print(m2);

    nml_mat_free(m1);
    nml_mat_free(m2);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {

  double v1[9] = {
    0.0, 1.0, 2.0,
    1.0, 2.0, 1.0,
    2.0, 7.0, 8.0
  };

  nml_mat *m1 = nml_mat_from(3, 3, 9, v1);
  printf("\nm1=\n");
  nml_mat_print(m1);
  nml_mat *refm1 = nml_mat_ref(m1);
  printf("\nrefm1=\n");
  nml_mat_print(refm1);

  nml_mat_free(m1);
  nml_mat_free(refm1);
  return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {

  double v1[9] = {
    0.0, 1.0, 2.0,
    1.0, 2.0, 1.0,
    2.0, 7.0, 8.0
  };

  nml_mat *m1 = nml_mat_from(3, 3, 9, v1);
  printf("\nm1=\n");
  nml_mat_print(m1);
  nml_mat *rrefm1 = nml_mat_rref(m1);
  printf("\nrrefm1=\n");
  nml_mat_print(rrefm1);

  nml_mat_free(m1);
  nml_mat_free(rrefm1);
  return 0;
}

typedef struct nml_mat_lup_s {
  nml_mat *L;
  nml_mat *U;
  nml_mat *P;
  unsigned int num_permutations;
} nml_mat_lup;

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *m1 = nml_mat_sqr_rnd(4, 0.0, 10.0);
    printf("m1=\n");
    nml_mat_print(m1);

    nml_mat_lup *m1_lup = nml_mat_lup_solve(m1);
    printf("L, U, P:\n");
    nml_mat_lup_print(m1_lup);

    nml_mat_free(m1);
    nml_mat_lup_free(m1_lup);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
  printf("\nInverse of a matrix:\n");
  double m_v[16] = {
    2.0, 7.0, 6.0, 1.0,
    9.0, 5.0, 0.0, 2.0,
    4.0, 3.0, 8.0, 3.0,
    3.0, 5.0, 1.0, 9.0
  };
  nml_mat *m = nml_mat_from(4,4,16, m_v);
  nml_mat_lup *lup = nml_mat_lup_solve(m);
  nml_mat* minv = nml_mat_inv(lup);
  nml_mat *mdotminv = nml_mat_dot(m, minv);

  printf("m=");
  nml_mat_print(m);
  printf("minv=");
  nml_mat_print(minv);
  printf("(%%e) m * minv=");
  nml_mat_printf(mdotminv, "%e\t");
  printf("(%%f) m * minv=");
  nml_mat_printf(mdotminv, "%f\t");

  nml_mat_free(m);
  nml_mat_free(minv);
  nml_mat_free(mdotminv);
  return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *m1 = nml_mat_sqr_rnd(4, 0.0, 10.0);
    nml_mat_lup *m1_lup = nml_mat_lup_solve(m1);
    
    printf("m1=\n");
    nml_mat_print(m1);
    printf("determinant=%lf\n", nml_mat_det(m1_lup));

    nml_mat_free(m1);
    nml_mat_lup_free(m1_lup);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    
    FILE *input = fopen("./examples/data/matrix9_lower_triangular.data", "r");

    nml_mat *A = nml_mat_fromfilef(input);
    nml_mat *B = nml_mat_fromfilef(input);
    nml_mat *x = nml_ls_solvefwd(A, B);

    nml_mat_print(A);
    nml_mat_print(B);
    nml_mat_print(x);

    nml_mat_free(A);
    nml_mat_free(B);
    nml_mat_free(x);

    fclose(input);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    
    FILE *input = fopen("./examples/data/matrix10_upper_triangular.data", "r");

    nml_mat *A = nml_mat_fromfilef(input);
    nml_mat *B = nml_mat_fromfilef(input);
    nml_mat *x = nml_ls_solvebck(A, B);

    nml_mat_print(A);
    nml_mat_print(B);
    nml_mat_print(x);

    nml_mat_free(A);
    nml_mat_free(B);
    nml_mat_free(x);

    fclose(input);

    return 0;
}

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

#include "lib/nml.h"

int main(int argc, char *argv[]) {
    nml_mat *A = nml_mat_sqr_rnd(4, 1.0, 10.0);
    nml_mat *B = nml_mat_rnd(4, 1, 1.0, 10.0);
    nml_mat_lup *LUP = nml_mat_lup_solve(A);

    nml_mat *x = nml_ls_solve(LUP, B);
    nml_mat_print(x);

    nml_mat_free(A);
    nml_mat_free(B);
    nml_mat_free(x);
    nml_mat_lup_free(LUP);
}